ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge pull request #1841 from phankl/stable

Browse Source 
Mesoscopic potential for carbon nanotubes
This commit is contained in:
Axel Kohlmeyer
2020-01-15 16:31:40 -05:00
committed by GitHub
parent 1ed50bcdde cb20d73a82
commit 3efc5b0b64
15 changed files with 2010084 additions and 1 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
doc/src/Commands_pair.rst
View File

@ -174,6 +174,7 @@ OPT.
* :doc:`meam/c <pair_meamc>`
* :doc:`meam/spline (o) <pair_meam_spline>`
* :doc:`meam/sw/spline <pair_meam_sw_spline>`
* :doc:`mesocnt <pair_mesocnt>`
* :doc:`mgpt <pair_mgpt>`
* :doc:`mie/cut (g) <pair_mie>`
* :doc:`mm3/switch3/coulgauss/long <pair_mm3_switch3_coulgauss_long>`
@ -257,4 +258,3 @@ OPT.
* :doc:`yukawa/colloid (go) <pair_yukawa_colloid>`
* :doc:`zbl (gko) <pair_zbl>`
*
*

161
doc/src/pair_mesocnt.rst Normal file
View File

@ -0,0 +1,161 @@
.. index:: pair\_style mesocnt
pair\_style mesocnt command
===========================
Syntax
""""""
.. parsed-literal::
pair_style mesocnt
Examples
""""""""
.. parsed-literal::
pair_style mesocnt
pair_coeff \* \* 10_10.cnt
Description
"""""""""""
Style *mesocnt* implements a mesoscopic potential
for the interaction of carbon nanotubes (CNTs). In this potential,
CNTs are modelled as chains of cylindrical segments in which
each infinitesimal surface element interacts with all other
CNT surface elements with the Lennard-Jones (LJ) term adopted from
the :doc:`airebo <pair_airebo>` style. The interaction energy
is then computed by integrating over the surfaces of all interacting
CNTs.
The potential is based on interactions between one cylindrical
segment and infinitely or semi-infinitely long CNTs as described
in :ref:`(Volkov1) <Volkov1>`. Chains of segments are
converted to these (semi-)infinite CNTs bases on an approximate
chain approach outlined in :ref:`(Volkov2) <Volkov2>`.
This allows to simplify the computation of the interactions
significantly and reduces the computational times to the
same order of magnitude as for regular bead spring models
where beads interact with the standard :doc:`pair_lj/cut <pair_lj>`
potential.
In LAMMPS, cylindrical segments are represented by bonds. Each
segment is defined by its two end points ("nodes") which correspond
to atoms in LAMMPS. For the exact functional form of the potential
and implementation details, the reader is referred to the
original papers :ref:`(Volkov1) <Volkov1>` and
:ref:`(Volkov2) <Volkov2>`.
The potential requires tabulated data provided in a single ASCII
text file specified in the :doc:`pair\_coeff <pair_coeff>` command.
The first line of the file provides a time stamp and
general information. The second line lists four integers giving
the number of data points provided in the subsequent four
data tables. The third line lists four floating point numbers:
the CNT radius R, the LJ parameter sigma and two numerical
parameters delta1 and delta2. These four parameters are given
in Angstroms. This is followed by four data tables each separated
by a single empty line. The first two tables have two columns
and list the parameters uInfParallel and Gamma respectively.
The last two tables have three columns giving data on a quadratic
array and list the parameters Phi and uSemiParallel respectively.
uInfParallel and uSemiParallel are given in eV/Angstrom, Phi is
given in eV and Gamma is unitless.
Potential files for CNTs can be readily generated using the freely
available code provided on
.. parsed-literal::
https://github.com/phankl/cntpot
Using the same approach, it should also be possible to
generate potential files for other 1D systems such as
boron nitride nanotubes.
.. note::
LAMMPS comes with one *mesocnt* style potential file
where the default number of data points per table is 1001.
This is sufficient for NVT simulations. For proper energy
conservation, we recommend using a potential file where
the resolution for Phi is at least 2001 data points.
.. note::
The *mesocnt* style requires CNTs to be represented
as a chain of atoms connected by bonds. Atoms need
to be numbered consecutively within one chain.
Atoms belonging to different CNTs need to be assigned
different molecule IDs.
A full summary of the method and LAMMPS implementation details
is expected to soon become available in Computer Physics
Communications.
----------
**Mixing, shift, table, tail correction, restart, rRESPA info**\ :
This pair style does not support mixing.
This pair style does not support the :doc:`pair\_modify <pair_modify>`
shift, table, and tail options.
The *mesocnt* pair style do not write their information to :doc:`binary restart files <restart>`,
since it is stored in tabulated potential files.
Thus, you need to re-specify the pair\_style and pair\_coeff commands in
an input script that reads a restart file.
This pair style can only be used via the *pair* keyword of the
:doc:`run\_style respa <run_style>` command. They do not support the
*inner*\ , *middle*\ , *outer* keywords.
----------
Restrictions
""""""""""""
This style is part of the USER-MISC package. It is only
enabled if LAMMPS was built with that package. See the :doc:`Build package <Build_package>` doc page for more info.
This pair potential requires the :doc:`newton <newton>` setting to be
"on" for pair interactions.
Related commands
""""""""""""""""
:doc:`pair\_coeff <pair_coeff>`
**Default:** none
----------
.. _Volkov1:
**(Volkov1)** Volkov and Zhigilei, J Phys Chem C, 114, 5513 (2010).
.. _Volkov2:
**(Volkov2)** Volkov, Simov and Zhigilei, APS Meeting Abstracts,
Q31.013 (2008).
.. _lws: http://lammps.sandia.gov
.. _ld: Manual.html
.. _lc: Commands_all.html

1
doc/src/pair_style.rst
View File

@ -244,6 +244,7 @@ accelerated styles exist.
* :doc:`meam/c <pair_meamc>` - modified embedded atom method (MEAM) in C
* :doc:`meam/spline <pair_meam_spline>` - splined version of MEAM
* :doc:`meam/sw/spline <pair_meam_sw_spline>` - splined version of MEAM with a Stillinger-Weber term
* :doc:`mesocnt <pair_mesocnt>` - mesoscale model for (carbon) nanotubes
* :doc:`mgpt <pair_mgpt>` - simplified model generalized pseudopotential theory (MGPT) potential
* :doc:`mie/cut <pair_mie>` - Mie potential
* :doc:`mm3/switch3/coulgauss/long <pair_mm3_switch3_coulgauss_long>` - smoothed MM3 vdW potential with Gaussian electrostatics

6
doc/utils/sphinx-config/false_positives.txt
View File

@ -1673,6 +1673,7 @@ Mersenne
Merz
meshless
meso
mesocnt
mesoparticle
mesoscale
mesoscopic
@ -2587,6 +2588,7 @@ Sikandar
Silbert
Silling
Sim
Simov
Simul
simulations
Sinnott
@ -2921,6 +2923,7 @@ UEF
ufm
Uhlenbeck
Ui
uInfParallel
uk
ul
ulb
@ -2961,6 +2964,7 @@ upto
Urbakh
Urbana
usec
uSemiParallel
userguide
username
usr
@ -3028,6 +3032,7 @@ VMDHOME
vn
Voigt
volfactor
Volkov
Volpe
volpress
volumetric
@ -3180,6 +3185,7 @@ zflag
Zhang
Zhen
zhi
Zhigilei
Zhou
Zhu
zi

1
examples/USER/misc/mesocnt/C_10_10.mesocnt Symbolic link
View File

@ -0,0 +1 @@
../../../../potentials/C_10_10.mesocnt

1521
examples/USER/misc/mesocnt/cnt.data Normal file
View File

File diff suppressed because it is too large Load Diff

38
examples/USER/misc/mesocnt/in.cnt Normal file
View File

@ -0,0 +1,38 @@
#Initialisation
units nano
dimension 3
boundary p p p
atom_style full
comm_modify cutoff 11.0
neighbor 7.80 bin
newton on
#Read data
read_data cnt.data
replicate 1 2 2
#Force field
bond_style harmonic
bond_coeff 1 268896.77 2.0
angle_style harmonic
angle_coeff 1 46562.17 180.0
pair_style mesocnt
pair_coeff * * C_10_10.mesocnt
#Output
thermo 1000
dump xyz all xyz 1000 cnt.xyz
#Simulation setup
timestep 1.0e-05
#Nose-Hoover thermostat
fix nvt all nvt temp 300 300 0.001
run 10000

126
examples/USER/misc/mesocnt/log.9Jan20.cnt.g++.1 Normal file
View File

@ -0,0 +1,126 @@
LAMMPS (09 Jan 2020)
#Initialisation
units nano
dimension 3
boundary p p p
atom_style full
comm_modify cutoff 11.0
neighbor 7.80 bin
newton on
#Read data
read_data cnt.data
orthogonal box = (0 0 0) to (600 600 60)
1 by 1 by 1 MPI processor grid
reading atoms ...
500 atoms
scanning bonds ...
1 = max bonds/atom
scanning angles ...
1 = max angles/atom
reading bonds ...
498 bonds
reading angles ...
496 angles
2 = max # of 1-2 neighbors
2 = max # of 1-3 neighbors
4 = max # of 1-4 neighbors
6 = max # of special neighbors
special bonds CPU = 0.000180006 secs
read_data CPU = 0.00125766 secs
replicate 1 2 2
orthogonal box = (0 0 0) to (600 1200 120)
1 by 1 by 1 MPI processor grid
2000 atoms
1992 bonds
1984 angles
2 = max # of 1-2 neighbors
2 = max # of 1-3 neighbors
4 = max # of 1-4 neighbors
6 = max # of special neighbors
special bonds CPU = 0.00054121 secs
replicate CPU = 0.000902414 secs
#Force field
bond_style harmonic
bond_coeff 1 268896.77 2.0
angle_style harmonic
angle_coeff 1 46562.17 180.0
pair_style mesocnt
pair_coeff * * 10_10.cnt
Reading potential file 10_10.cnt with DATE: 2020-01-13
#Output
thermo 1000
dump xyz all xyz 1000 cnt.xyz
#Simulation setup
timestep 1.0e-05
#Nose-Hoover thermostat
fix nvt all nvt temp 300 300 0.001
run 10000
WARNING: Using a manybody potential with bonds/angles/dihedrals and special_bond exclusions (src/pair.cpp:226)
Neighbor list info ...
update every 1 steps, delay 10 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 10.177
ghost atom cutoff = 11
binsize = 5.0885, bins = 118 236 24
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair mesocnt, perpetual
attributes: full, newton on
pair build: full/bin
stencil: full/bin/3d
bin: standard
Per MPI rank memory allocation (min/avg/max) = 11.21 | 11.21 | 11.21 Mbytes
Step Temp E_pair E_mol TotEng Press
0 0 -4632.0136 5.7939238e-17 -4632.0136 -0.00015032304
1000 298.82235 -19950.274 13208.089 5628.7024 -0.0056205182
2000 300.43933 -28320.212 11980.296 -3902.0877 -0.0045324757
3000 300.4263 -36049.855 11338.405 -12274.161 -0.0018833539
4000 299.13368 -43471.21 11926.882 -19160.553 -0.00043030866
5000 293.77858 -50083.893 12334.927 -25586.884 -0.0015653738
6000 296.4851 -56330.135 12325.63 -31730.376 -0.0012795986
7000 298.20879 -62120.359 12582.297 -37192.574 -0.0013845796
8000 299.45547 -67881.692 13058.926 -42425.669 -0.00021100885
9000 301.82622 -73333.698 13598.257 -47240.197 -0.0006009197
10000 307.16873 -78292.306 13818.929 -51756.96 -0.0005609903
Loop time of 4.0316 on 1 procs for 10000 steps with 2000 atoms
Performance: 2143.072 ns/day, 0.011 hours/ns, 2480.408 timesteps/s
99.9% CPU use with 1 MPI tasks x no OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 2.5955 | 2.5955 | 2.5955 | 0.0 | 64.38
Bond | 1.1516 | 1.1516 | 1.1516 | 0.0 | 28.57
Neigh | 0.001163 | 0.001163 | 0.001163 | 0.0 | 0.03
Comm | 0.0019577 | 0.0019577 | 0.0019577 | 0.0 | 0.05
Output | 0.020854 | 0.020854 | 0.020854 | 0.0 | 0.52
Modify | 0.21637 | 0.21637 | 0.21637 | 0.0 | 5.37
Other | | 0.04409 | | | 1.09
Nlocal: 2000 ave 2000 max 2000 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 0 ave 0 max 0 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 0 ave 0 max 0 min
Histogram: 1 0 0 0 0 0 0 0 0 0
FullNghs: 13320 ave 13320 max 13320 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 13320
Ave neighs/atom = 6.66
Ave special neighs/atom = 5.952
Neighbor list builds = 1
Dangerous builds = 0
Total wall time: 0:00:05

126
examples/USER/misc/mesocnt/log.9Jan20.cnt.g++.4 Normal file
View File

@ -0,0 +1,126 @@
LAMMPS (09 Jan 2020)
#Initialisation
units nano
dimension 3
boundary p p p
atom_style full
comm_modify cutoff 11.0
neighbor 7.80 bin
newton on
#Read data
read_data cnt.data
orthogonal box = (0 0 0) to (600 600 60)
2 by 2 by 1 MPI processor grid
reading atoms ...
500 atoms
scanning bonds ...
1 = max bonds/atom
scanning angles ...
1 = max angles/atom
reading bonds ...
498 bonds
reading angles ...
496 angles
2 = max # of 1-2 neighbors
2 = max # of 1-3 neighbors
4 = max # of 1-4 neighbors
6 = max # of special neighbors
special bonds CPU = 0.000354767 secs
read_data CPU = 0.00286365 secs
replicate 1 2 2
orthogonal box = (0 0 0) to (600 1200 120)
1 by 4 by 1 MPI processor grid
2000 atoms
1992 bonds
1984 angles
2 = max # of 1-2 neighbors
2 = max # of 1-3 neighbors
4 = max # of 1-4 neighbors
6 = max # of special neighbors
special bonds CPU = 0.00019598 secs
replicate CPU = 0.00055337 secs
#Force field
bond_style harmonic
bond_coeff 1 268896.77 2.0
angle_style harmonic
angle_coeff 1 46562.17 180.0
pair_style mesocnt
pair_coeff * * 10_10.cnt
Reading potential file 10_10.cnt with DATE: 2020-01-13
#Output
thermo 1000
dump xyz all xyz 1000 cnt.xyz
#Simulation setup
timestep 1.0e-05
#Nose-Hoover thermostat
fix nvt all nvt temp 300 300 0.001
run 10000
WARNING: Using a manybody potential with bonds/angles/dihedrals and special_bond exclusions (src/pair.cpp:226)
Neighbor list info ...
update every 1 steps, delay 10 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 10.177
ghost atom cutoff = 11
binsize = 5.0885, bins = 118 236 24
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair mesocnt, perpetual
attributes: full, newton on
pair build: full/bin
stencil: full/bin/3d
bin: standard
Per MPI rank memory allocation (min/avg/max) = 2.725 | 2.725 | 2.725 Mbytes
Step Temp E_pair E_mol TotEng Press
0 0 -4632.0136 5.7939238e-17 -4632.0136 -0.00015032304
1000 298.82235 -19950.274 13208.089 5628.7024 -0.0056205182
2000 300.43861 -28320.205 11980.287 -3902.1202 -0.0045324738
3000 300.41076 -36049.308 11339.149 -12273.513 -0.0018848513
4000 299.13326 -43471.424 11927.668 -19159.998 -0.00042845101
5000 293.78857 -50083.216 12333.969 -25586.752 -0.0015664633
6000 296.45482 -56329.621 12326.419 -31730.328 -0.0012773686
7000 298.19097 -62119.086 12581.4 -37192.937 -0.0013862831
8000 299.46424 -67880.989 13057.62 -42425.908 -0.00020874264
9000 301.80677 -73332.208 13597.237 -47240.532 -0.00060074773
10000 307.17104 -78292.912 13818.889 -51757.51 -0.00056148282
Loop time of 1.23665 on 4 procs for 10000 steps with 2000 atoms
Performance: 6986.607 ns/day, 0.003 hours/ns, 8086.351 timesteps/s
96.1% CPU use with 4 MPI tasks x no OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 0.66321 | 0.68439 | 0.71413 | 2.5 | 55.34
Bond | 0.28561 | 0.29434 | 0.30976 | 1.7 | 23.80
Neigh | 0.00043321 | 0.00043637 | 0.00043917 | 0.0 | 0.04
Comm | 0.026656 | 0.05346 | 0.097228 | 12.7 | 4.32
Output | 0.0070224 | 0.0073031 | 0.0081415 | 0.6 | 0.59
Modify | 0.12769 | 0.15394 | 0.18743 | 6.5 | 12.45
Other | | 0.04279 | | | 3.46
Nlocal: 500 ave 504 max 496 min
Histogram: 2 0 0 0 0 0 0 0 0 2
Nghost: 22 ave 24 max 20 min
Histogram: 2 0 0 0 0 0 0 0 0 2
Neighs: 0 ave 0 max 0 min
Histogram: 4 0 0 0 0 0 0 0 0 0
FullNghs: 3330 ave 3368 max 3292 min
Histogram: 2 0 0 0 0 0 0 0 0 2
Total # of neighbors = 13320
Ave neighs/atom = 6.66
Ave special neighs/atom = 5.952
Neighbor list builds = 1
Dangerous builds = 0
Total wall time: 0:00:02

2006011
potentials/C_10_10.mesocnt Normal file
View File

File diff suppressed because it is too large Load Diff

1
potentials/README
View File

@ -89,6 +89,7 @@ lcbop LCBOP long-range bond-order potential
meam modified EAM (MEAM) library and individual elements/alloys
meam.spline modified EAM (MEAM) spline potential
meam.sw.spline modified EAM (MEAM) Stillinger-Weber spline potential
mesocnt mesoscopic carbon nanotube (CNT) potential
mgpt model generalized pseudopotential theory (MGPT) potential
nb3b.harmonic nonbonded 3-body harmonic potential
poly polymorphic 3-body potential

1
src/USER-MISC/README
View File

@ -91,6 +91,7 @@ pair_style kolmogorov/crespi/full, Wengen Ouyang (Tel Aviv University), w.g.ouya
pair_style kolmogorov/crespi/z, Jaap Kroes (Radboud U), jaapkroes at gmail dot com, 28 Feb 17
pair_style meam/spline, Alexander Stukowski (LLNL), alex at stukowski.com, 1 Feb 12
pair_style meam/sw/spline, Robert Rudd (LLNL), robert.rudd at llnl.gov, 1 Oct 12
pair_style mesocnt, Philipp Kloza (U Cambridge), pak37 at cam.ac.uk, 15 Jan 20
pair_style morse/smooth/linear, Stefan Paquay (TU Eindhoven), stefanpaquay at gmail.com, 29 Feb 16
pair_style srp, Tim Sirk, tim.sirk at us.army.mil, 21 Nov 14
pair_style tersoff/table, Luca Ferraro, luca.ferraro@caspur.it, 1 Dec 11

1837
src/USER-MISC/pair_mesocnt.cpp Normal file
View File

File diff suppressed because it is too large Load Diff

170
src/USER-MISC/pair_mesocnt.h Normal file
View File

@ -0,0 +1,170 @@
#ifdef PAIR_CLASS
PairStyle(mesocnt, PairMesoCNT)
#else
#ifndef LMP_PAIR_MESOCNT_H
#define LMP_PAIR_MESOCNT_H
#include "pair.h"
namespace LAMMPS_NS {
class PairMesoCNT : public Pair {
public:
PairMesoCNT(class LAMMPS *);
~PairMesoCNT();
void compute(int, int);
void settings(int, char **);
void coeff(int, char **);
void init_style();
double init_one(int, int);
protected:
int uinf_points,gamma_points,phi_points,usemi_points;
int nlocal_size,reduced_neigh_size;
int *reduced_nlist,*numchainlist;
int **reduced_neighlist,**nchainlist,**endlist;
int ***chainlist;
double ang,ang_inv,eunit,funit;
double delta1,delta2;
double r,rsq,d,rc,rcsq,rc0,cutoff,cutoffsq;
double r_ang,rsq_ang,d_ang,rc_ang,rcsq_ang,cutoff_ang,cutoffsq_ang;
double sig,sig_ang,comega,ctheta;
double hstart_uinf,hstart_gamma,
hstart_phi,psistart_phi,hstart_usemi,xistart_usemi;
double delh_uinf,delh_gamma,delh_phi,delpsi_phi,delh_usemi,delxi_usemi;
double p1[3],p2[3],p[3],m[3];
double *param,*w,*wnode;
double **dq_w;
double ***q1_dq_w,***q2_dq_w;
double *uinf_data,*gamma_data,**phi_data,**usemi_data;
double **uinf_coeff,**gamma_coeff,****phi_coeff,****usemi_coeff;
double **flocal,**fglobal,**basis;
char *file;
void allocate();
void bond_neigh();
void neigh_common(int, int, int &, int *);
void chain_split(int *, int, int &, int **, int *, int *);
void sort(int *, int);
void read_file();
void read_data(FILE *, double *, double &, double &, int);
void read_data(FILE *, double **, double &, double &,
double &, double &, int);
void spline_coeff(double *, double **, double, int);
void spline_coeff(double **, double ****, double, double, int);
double spline(double, double, double, double **, int);
double dspline(double, double, double, double **, int);
double spline(double, double, double, double, double, double,
double ****, int);
double dxspline(double, double, double, double, double, double,
double ****, int);
double dyspline(double, double, double, double, double, double,
double ****, int);
void geometry(const double *, const double *, const double *,
const double *, const double *,
double *, double *, double *, double **);
void weight(const double *, const double *, const double *,
const double *, double &, double *, double *,
double *, double *);
void finf(const double *, double &, double **);
void fsemi(const double *, double &, double &, double **);
// inlined functions for efficiency
inline double heaviside(double x) {
if (x > 0) return 1.0;
else return 0.0;
}
inline double s(double x) {
return heaviside(-x) + heaviside(x)*heaviside(1-x)*(1 - x*x*(3 - 2*x));
}
inline double ds(double x) {
return 6 * heaviside(x) * heaviside(1-x) * x * (x-1);
}
inline double s5(double x) {
double x2 = x * x;
return heaviside(-x)
+ heaviside(x)*heaviside(1-x)*(1 - x2*x*(6*x2 - 15*x + 10));
}
inline double ds5(double x) {
double x2 = x * x;
return -30 * heaviside(x) * heaviside(1-x) * x2 * (x2 - 2*x + 1);
}
};
}
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Incorrect args for pair coefficients
Self-explanatory. Check the input script or data file.
E: Pair style mesocnt does not support lj units
Self-explanatory. Specify different unit system using the units
command.
E: Pair style mesocnt requires atom IDs
Self-explanatory. Turn on atom IDs using the atom_modify command.
E: Pair style mesocnt requires newton pair on
Self-explanatory. Turn on Newton's third law with the newton command.
E: Cannot open mesocnt file %s
The specified mesocnt potential file cannot be opened. Check that the
path and name are correct.
E: Premature end of file in pair table %s
The specified mesocnt potential file is shorter than specified. Check
if the correct file is being used and the right number of data points
was specified in the pair_style
W: %d of %d lines were incomplete or could not be parsed completely
in pair table %s
A number of lines in the specified mesocnt potential file is incomplete
or in the wrong format. Check the file for errors and missing data.
W: %d spacings in the first column were different from the first spacing
in the pair table %s
The spacings between x coordinates in the first column of the specified
mesocnt potential file vary throughout the file. Use a potential file
with higher precision.
W: %d spacings in second column were different from first
spacing in pair table %s
The spacings between y coordinates in the second column of the specified
mesocnt potential file vary throughout the file. Use a potential file
with higher precision.
*/

83
src/math_extra.h
View File

@ -43,6 +43,9 @@ namespace MathExtra {
// 3x3 matrix operations
inline void zeromat3(double m[3][3]);
inline void zeromat3(double **m);
inline void col2mat(const double *ex, const double *ey, const double *ez,
double m[3][3]);
inline double det3(const double mat[3][3]);
@ -52,6 +55,12 @@ namespace MathExtra {
double ans[3][3]);
inline void plus3(const double m[3][3], const double m2[3][3],
double ans[3][3]);
inline void plus3(const double m[3][3], double **m2, double **ans);
inline void minus3(const double m[3][3], const double m2[3][3],
double ans[3][3]);
inline void minus3(double **m, const double m2[3][3],
double ans[3][3]);
inline void times3(const double m[3][3], const double m2[3][3],
double ans[3][3]);
inline void transpose_times3(const double m[3][3], const double m2[3][3],
@ -71,6 +80,8 @@ namespace MathExtra {
double ans[3][3]);
inline void vecmat(const double *v, const double m[3][3], double *ans);
inline void scalar_times3(const double f, double m[3][3]);
inline void outer3(const double *v1, const double *v2,
double ans[3][3]);
void write3(const double mat[3][3]);
int mldivide3(const double mat[3][3], const double *vec, double *ans);
@ -707,4 +718,76 @@ inline void MathExtra::rotation_generator_z(const double m[3][3],
ans[2][2] = 0;
}
// set matrix to zero
inline void MathExtra::zeromat3(double m[3][3])
{
m[0][0] = m[0][1] = m[0][2] = 0.0;
m[1][0] = m[1][1] = m[1][2] = 0.0;
m[2][0] = m[2][1] = m[2][2] = 0.0;
}
inline void MathExtra::zeromat3(double **m)
{
m[0][0] = m[0][1] = m[0][2] = 0.0;
m[1][0] = m[1][1] = m[1][2] = 0.0;
m[2][0] = m[2][1] = m[2][2] = 0.0;
}
// add two matrices
inline void MathExtra::plus3(const double m[3][3], double **m2,
double **ans)
{
ans[0][0] = m[0][0]+m2[0][0];
ans[0][1] = m[0][1]+m2[0][1];
ans[0][2] = m[0][2]+m2[0][2];
ans[1][0] = m[1][0]+m2[1][0];
ans[1][1] = m[1][1]+m2[1][1];
ans[1][2] = m[1][2]+m2[1][2];
ans[2][0] = m[2][0]+m2[2][0];
ans[2][1] = m[2][1]+m2[2][1];
ans[2][2] = m[2][2]+m2[2][2];
}
// subtract two matrices
inline void MathExtra::minus3(const double m[3][3], const double m2[3][3],
double ans[3][3])
{
ans[0][0] = m[0][0]-m2[0][0];
ans[0][1] = m[0][1]-m2[0][1];
ans[0][2] = m[0][2]-m2[0][2];
ans[1][0] = m[1][0]-m2[1][0];
ans[1][1] = m[1][1]-m2[1][1];
ans[1][2] = m[1][2]-m2[1][2];
ans[2][0] = m[2][0]-m2[2][0];
ans[2][1] = m[2][1]-m2[2][1];
ans[2][2] = m[2][2]-m2[2][2];
}
inline void MathExtra::minus3(double **m, const double m2[3][3],
double ans[3][3])
{
ans[0][0] = m[0][0]-m2[0][0];
ans[0][1] = m[0][1]-m2[0][1];
ans[0][2] = m[0][2]-m2[0][2];
ans[1][0] = m[1][0]-m2[1][0];
ans[1][1] = m[1][1]-m2[1][1];
ans[1][2] = m[1][2]-m2[1][2];
ans[2][0] = m[2][0]-m2[2][0];
ans[2][1] = m[2][1]-m2[2][1];
ans[2][2] = m[2][2]-m2[2][2];
}
// compute outer product of two vectors
inline void MathExtra::outer3(const double *v1, const double *v2,
double ans[3][3])
{
ans[0][0] = v1[0]*v2[0]; ans[0][1] = v1[0]*v2[1]; ans[0][2] = v1[0]*v2[2];
ans[1][0] = v1[1]*v2[0]; ans[1][1] = v1[1]*v2[1]; ans[1][2] = v1[1]*v2[2];
ans[2][0] = v1[2]*v2[0]; ans[2][1] = v1[2]*v2[1]; ans[2][2] = v1[2]*v2[2];
}
#endif