diff --git a/doc/src/compute_born_matrix.rst b/doc/src/compute_born_matrix.rst new file mode 100644 index 0000000000..2a9fd46e1f --- /dev/null +++ b/doc/src/compute_born_matrix.rst @@ -0,0 +1,136 @@ +.. index:: compute born/matrix + +compute born/matrix command +=========================== + +Syntax +"""""" + +.. parsed-literal:: + + compute ID group-ID born/matrix + +* ID, group-ID are documented in :doc:`compute ` command +* born/matrix = style name of this compute command + +Examples +"""""""" + +.. code-block:: LAMMPS + + compute 1 all born/matrix + +Description +""""""""""" + +Define a compute that calculates +:math:`\frac{\partial{}^2U}{\partial\varepsilon_{i}\partial\varepsilon_{j}}` the +second derivatives of the potential energy :math:`U` with regard to strain +tensor :math:`\varepsilon` elements. These values are related to: + +.. math:: + + C^{B}_{i,j}=\frac{1}{V}\frac{\partial{}^2U}{\partial{}\varepsilon_{i}\partial\varepsilon_{j}} + +also called the Born term of elastic constants in the stress-stress fluctuation +formalism. This quantity can be used to compute the elastic constant tensor. +Using the symmetric Voigt notation, the elastic constant tensor can be written +as a 6x6 symmetric matrix: + +.. math:: + + C_{i,j} = \langle{}C^{B}_{i,j}\rangle + + \frac{V}{k_{B}T}\left(\langle\sigma_{i}\sigma_{j}\rangle\right. + \left.- \langle\sigma_{i}\rangle\langle\sigma_{j}\rangle\right) + + \frac{Nk_{B}T}{V} + \left(\delta_{i,j}+(\delta_{1,i}+\delta_{2,i}+\delta_{3,i})\right. + \left.*(\delta_{1,j}+\delta_{2,j}+\delta_{3,j})\right) + +In the above expression, :math:`\sigma` stands for the virial stress +tensor, :math:`\delta` is the Kronecker delta and the usual notation apply for +the number of particle, the temperature and volume respectively :math:`N`, +:math:`T` and :math:`V`. :math:`k_{B}` is the Boltzmann constant. + +The Born term is a symmetric 6x6 matrix by construction and as such can be +expressed as 21 independent terms. The terms are ordered corresponding to the +following matrix element: + +.. math:: + + \begin{matrix} + C_{1} & C_{7} & C_{8} & C_{9} & C_{10} & C_{11} \\ + C_{7} & C_{2} & C_{12} & C_{13} & C_{14} & C_{15} \\ + \vdots & C_{12} & C_{3} & C_{16} & C_{17} & C_{18} \\ + \vdots & C_{13} & C_{16} & C_{4} & C_{19} & C_{20} \\ + \vdots & \vdots & \vdots & C_{19} & C_{5} & C_{21} \\ + \vdots & \vdots & \vdots & \vdots & C_{21} & C_{6} + \end{matrix} + +in this matrix the indices of :math:`C_{k}` value are the corresponding index +:math:`k` in the compute output. Each term comes from the sum of every +interactions derivatives in the system as explained in :ref:`(VanWorkum) +` or :ref:`(Voyiatzis) `. + +The output can be accessed using usual Lammps routines: + +.. code-block:: LAMMPS + + compute 1 all born/matrix + compute 2 all pressure NULL virial + variable S1 equal -c_2[1] + variable S2 equal -c_2[2] + variable S3 equal -c_2[3] + variable S4 equal -c_2[4] + variable S5 equal -c_2[5] + variable S6 equal -c_2[6] + fix 1 all ave/time 1 1 1 v_S1 v_S2 v_S3 v_S4 v_S5 v_S6 c_1[*] file born.out + +In this example, the file *born.out* will contain the information needed to +compute the first and second terms of the elastic constant matrix in a post +processing procedure. The other required quantities can be accessed using any +other *LAMMPS* usual method. + +NOTE: In the above :math:`C_{i,j}` computation, the term involving the virial +stress tensor :math:`\sigma` is the covariance between each elements. In a +solid the virial stress can have large variations between timesteps and average +values can be slow to converge. This term is better computed using +instantaneous values. + +**Output info:** + +This compute calculates a global array with the number of rows=21. +The values are ordered as explained above. These values can be used +by any command that uses a global values from a compute as input. See +the :doc:`Howto output ` doc page for an overview of +LAMMPS output options. + +The array values calculated by this compute are all "extensive". + +Restrictions +"""""""""""" + +This compute is part of the EXTRA-COMPUTE package. It is only enabled if +LAMMPS was built with that package. See the :doc:`Build package ` page for more info. + +The Born term can be decomposed as a product of two terms. The first one +is a general term which depends on the configuration. The second one is +specific to every interaction composing your force field (non-bonded, +bonds, angle...). Currently not all interaction implement the *born_matrix* +method giving first and second order derivatives and a warning will +be raised if you try to use this compute with such interactions. The returned +values of this force field component is currently zero. + +Default +""""""" + +none + +---------- + +.. _VanWorkum: + +**(Van Workum)** K. Van Workum et al., J. Chem. Phys. 125 144506 (2006) + +.. _Voyiatzis: + +**(Voyiatzis)** E. Voyiatzis, Computer Physics Communications 184(2013)27-33 diff --git a/examples/cij_nvt/in.ljfcc b/examples/cij_nvt/in.ljfcc new file mode 100644 index 0000000000..71e541df24 --- /dev/null +++ b/examples/cij_nvt/in.ljfcc @@ -0,0 +1,36 @@ +# Numerical difference calculation +# of Born matrix + +# adjustable parameters + +variable nsteps index 50 # length of run +variable nthermo index 10 # thermo output interval +variable nlat index 3 # size of box +variable delta index 1.0e-6 # strain size +variable temp index 0.3 # temperature + +units lj +atom_style atomic + +lattice fcc 0.8442 +region box block 0 ${nlat} 0 ${nlat} 0 ${nlat} +create_box 1 box +create_atoms 1 box +mass 1 1.0 + +velocity all create ${temp} 87287 loop geom + +pair_style lj/cut 2.5 +pair_coeff 1 1 1.0 1.0 + +neighbor 0.0 bin +neigh_modify every 1 delay 0 check no + +fix 1 all nve + +compute born all born/matrix +thermo ${nthermo} +thermo_style custom step temp pe press c_born[*] + +run ${nsteps} + diff --git a/examples/cij_nvt/in.ljfcc_num b/examples/cij_nvt/in.ljfcc_num new file mode 100644 index 0000000000..9f590cc0fd --- /dev/null +++ b/examples/cij_nvt/in.ljfcc_num @@ -0,0 +1,37 @@ +# Numerical difference calculation +# of Born matrix + +# adjustable parameters + +variable nsteps index 50 # length of run +variable nthermo index 10 # thermo output interval +variable nlat index 3 # size of box +variable delta index 1.0e-6 # strain size +variable temp index 0.3 # temperature + +units lj +atom_style atomic + +lattice fcc 0.8442 +region box block 0 ${nlat} 0 ${nlat} 0 ${nlat} +create_box 1 box +create_atoms 1 box +mass 1 1.0 + +velocity all create ${temp} 87287 loop geom + +pair_style lj/cut 2.5 +pair_coeff 1 1 1.0 1.0 + +neighbor 0.0 bin +neigh_modify every 1 delay 0 check no + +fix 1 all nve + +compute myvirial all pressure NULL virial +compute born all born/matrix numdiff ${delta} myvirial +thermo ${nthermo} +thermo_style custom step temp pe press c_born[*] + +run ${nsteps} + diff --git a/examples/cij_nvt/log.29Jan2022.born.g++.1 b/examples/cij_nvt/log.29Jan2022.born.g++.1 new file mode 100644 index 0000000000..b40a3e8676 --- /dev/null +++ b/examples/cij_nvt/log.29Jan2022.born.g++.1 @@ -0,0 +1,102 @@ +LAMMPS (7 Jan 2022) +# Numerical difference calculation +# of Born matrix + +# adjustable parameters + +variable nsteps index 50 # length of run +variable nthermo index 10 # thermo output interval +variable nlat index 3 # size of box +variable delta index 1.0e-6 # strain size +variable temp index 0.3 # temperature + +units lj +atom_style atomic + +lattice fcc 0.8442 +Lattice spacing in x,y,z = 1.6795962 1.6795962 1.6795962 +region box block 0 ${nlat} 0 ${nlat} 0 ${nlat} +region box block 0 3 0 ${nlat} 0 ${nlat} +region box block 0 3 0 3 0 ${nlat} +region box block 0 3 0 3 0 3 +create_box 1 box +Created orthogonal box = (0 0 0) to (5.0387886 5.0387886 5.0387886) + 1 by 1 by 1 MPI processor grid +create_atoms 1 box +Created 108 atoms + using lattice units in orthogonal box = (0 0 0) to (5.0387886 5.0387886 5.0387886) + create_atoms CPU = 0.000 seconds +mass 1 1.0 + +velocity all create ${temp} 87287 loop geom +velocity all create 0.3 87287 loop geom + +pair_style lj/cut 2.5 +pair_coeff 1 1 1.0 1.0 + +neighbor 0.0 bin +neigh_modify every 1 delay 0 check no + +fix 1 all nve + +compute born all born/matrix +thermo ${nthermo} +thermo 10 +thermo_style custom step temp pe press c_born[*] + +run ${nsteps} +run 50 + generated 0 of 0 mixed pair_coeff terms from geometric mixing rule +Neighbor list info ... + update every 1 steps, delay 0 steps, check no + max neighbors/atom: 2000, page size: 100000 + master list distance cutoff = 2.5 + ghost atom cutoff = 2.5 + binsize = 1.25, bins = 5 5 5 + 2 neighbor lists, perpetual/occasional/extra = 1 1 0 + (1) pair lj/cut, perpetual + attributes: half, newton on + pair build: half/bin/atomonly/newton + stencil: half/bin/3d + bin: standard + (2) compute born/matrix, occasional, copy from (1) + attributes: half, newton on + pair build: copy + stencil: none + bin: none +Per MPI rank memory allocation (min/avg/max) = 2.574 | 2.574 | 2.574 Mbytes +Step Temp PotEng Press c_born[1] c_born[2] c_born[3] c_born[4] c_born[5] c_born[6] c_born[7] c_born[8] c_born[9] c_born[10] c_born[11] c_born[12] c_born[13] c_born[14] c_born[15] c_born[16] c_born[17] c_born[18] c_born[19] c_born[20] c_born[21] + 0 0.3 -6.7733681 -5.9844023 411.5161 411.5161 411.5161 595.57522 595.57522 595.57522 595.57522 595.57522 -2.7200464e-15 -6.1062266e-15 7.7715612e-16 595.57522 -1.1518564e-14 -3.3306691e-15 -5.3013149e-15 -1.0630385e-14 -6.800116e-15 -2.6367797e-15 -3.4139358e-15 -3.7747583e-15 -2.7200464e-15 + 10 0.25283594 -6.703292 -5.5569889 891.93989 928.15437 944.55768 872.69866 852.24088 816.0675 816.0675 852.24088 10.111334 3.6921122 -6.3839362 872.69866 -4.6709643 4.4335317 -9.7018557 -4.5201548 3.7624526 11.200776 11.200776 4.4335317 10.111334 + 20 0.13298345 -6.5218396 -4.5683635 1958.2215 2130.5784 2172.9636 1516.3032 1426.6019 1323.0834 1323.0834 1426.6019 39.867158 11.406477 -25.845853 1516.3032 3.3483658 27.179119 -36.809191 10.147948 7.0018777 41.995438 41.995438 27.179119 39.867158 + 30 0.089302796 -6.4604227 -4.1807379 2401.2346 2625.1676 2638.9312 1724.973 1625.4639 1602.8817 1602.8817 1625.4639 55.594925 -101.03316 -41.689385 1724.973 35.280816 26.589804 -40.046683 63.973353 -78.086653 41.559746 41.559746 26.589804 55.594925 + 40 0.12149354 -6.5076059 -4.4000835 2234.612 2345.8821 2337.9552 1513.845 1486.9674 1564.2263 1564.2263 1486.9674 56.775913 -199.0226 -45.302207 1513.845 27.333643 -9.5839331 -16.768077 120.5015 -141.06684 17.802047 17.802047 -9.5839331 56.775913 + 50 0.13053032 -6.5204024 -4.5064456 2097.1985 2270.3252 2180.7348 1458.5923 1396.6336 1507.4402 1507.4402 1396.6336 37.603097 -213.85355 22.82129 1458.5923 47.47204 -20.546302 54.328025 195.36071 -169.33495 -24.624815 -24.624815 -20.546302 37.603097 +Loop time of 0.00790344 on 1 procs for 50 steps with 108 atoms + +Performance: 2732988.192 tau/day, 6326.362 timesteps/s +98.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 | 0.0014593 | 0.0014593 | 0.0014593 | 0.0 | 18.46 +Neigh | 0.0047863 | 0.0047863 | 0.0047863 | 0.0 | 60.56 +Comm | 0.00057975 | 0.00057975 | 0.00057975 | 0.0 | 7.34 +Output | 0.00096954 | 0.00096954 | 0.00096954 | 0.0 | 12.27 +Modify | 5.3432e-05 | 5.3432e-05 | 5.3432e-05 | 0.0 | 0.68 +Other | | 5.508e-05 | | | 0.70 + +Nlocal: 108 ave 108 max 108 min +Histogram: 1 0 0 0 0 0 0 0 0 0 +Nghost: 700 ave 700 max 700 min +Histogram: 1 0 0 0 0 0 0 0 0 0 +Neighs: 2909 ave 2909 max 2909 min +Histogram: 1 0 0 0 0 0 0 0 0 0 + +Total # of neighbors = 2909 +Ave neighs/atom = 26.935185 +Neighbor list builds = 50 +Dangerous builds not checked + +Total wall time: 0:00:00 diff --git a/examples/cij_nvt/log.29Jan2022.born.g++.4 b/examples/cij_nvt/log.29Jan2022.born.g++.4 new file mode 100644 index 0000000000..e5e2c207d6 --- /dev/null +++ b/examples/cij_nvt/log.29Jan2022.born.g++.4 @@ -0,0 +1,102 @@ +LAMMPS (7 Jan 2022) +# Numerical difference calculation +# of Born matrix + +# adjustable parameters + +variable nsteps index 50 # length of run +variable nthermo index 10 # thermo output interval +variable nlat index 3 # size of box +variable delta index 1.0e-6 # strain size +variable temp index 0.3 # temperature + +units lj +atom_style atomic + +lattice fcc 0.8442 +Lattice spacing in x,y,z = 1.6795962 1.6795962 1.6795962 +region box block 0 ${nlat} 0 ${nlat} 0 ${nlat} +region box block 0 3 0 ${nlat} 0 ${nlat} +region box block 0 3 0 3 0 ${nlat} +region box block 0 3 0 3 0 3 +create_box 1 box +Created orthogonal box = (0 0 0) to (5.0387886 5.0387886 5.0387886) + 1 by 2 by 2 MPI processor grid +create_atoms 1 box +Created 108 atoms + using lattice units in orthogonal box = (0 0 0) to (5.0387886 5.0387886 5.0387886) + create_atoms CPU = 0.000 seconds +mass 1 1.0 + +velocity all create ${temp} 87287 loop geom +velocity all create 0.3 87287 loop geom + +pair_style lj/cut 2.5 +pair_coeff 1 1 1.0 1.0 + +neighbor 0.0 bin +neigh_modify every 1 delay 0 check no + +fix 1 all nve + +compute born all born/matrix +thermo ${nthermo} +thermo 10 +thermo_style custom step temp pe press c_born[*] + +run ${nsteps} +run 50 + generated 0 of 0 mixed pair_coeff terms from geometric mixing rule +Neighbor list info ... + update every 1 steps, delay 0 steps, check no + max neighbors/atom: 2000, page size: 100000 + master list distance cutoff = 2.5 + ghost atom cutoff = 2.5 + binsize = 1.25, bins = 5 5 5 + 2 neighbor lists, perpetual/occasional/extra = 1 1 0 + (1) pair lj/cut, perpetual + attributes: half, newton on + pair build: half/bin/atomonly/newton + stencil: half/bin/3d + bin: standard + (2) compute born/matrix, occasional, copy from (1) + attributes: half, newton on + pair build: copy + stencil: none + bin: none +Per MPI rank memory allocation (min/avg/max) = 2.561 | 2.561 | 2.561 Mbytes +Step Temp PotEng Press c_born[1] c_born[2] c_born[3] c_born[4] c_born[5] c_born[6] c_born[7] c_born[8] c_born[9] c_born[10] c_born[11] c_born[12] c_born[13] c_born[14] c_born[15] c_born[16] c_born[17] c_born[18] c_born[19] c_born[20] c_born[21] + 0 0.3 -6.7733681 -5.9844023 411.5161 411.5161 411.5161 595.57522 595.57522 595.57522 595.57522 595.57522 -2.9143354e-15 -1.9984014e-15 1.8873791e-15 595.57522 -9.7560848e-15 -3.3029135e-15 -9.6311847e-15 -1.2281842e-14 -5.3568261e-15 -2.5257574e-15 -3.2196468e-15 -3.5804693e-15 -2.553513e-15 + 10 0.25283594 -6.703292 -5.5569889 891.93989 928.15437 944.55768 872.69866 852.24088 816.0675 816.0675 852.24088 10.111334 3.6921122 -6.3839362 872.69866 -4.6709643 4.4335317 -9.7018557 -4.5201548 3.7624526 11.200776 11.200776 4.4335317 10.111334 + 20 0.13298345 -6.5218396 -4.5683635 1958.2215 2130.5784 2172.9636 1516.3032 1426.6019 1323.0834 1323.0834 1426.6019 39.867158 11.406477 -25.845853 1516.3032 3.3483658 27.179119 -36.809191 10.147948 7.0018777 41.995438 41.995438 27.179119 39.867158 + 30 0.089302796 -6.4604227 -4.1807379 2401.2346 2625.1676 2638.9312 1724.973 1625.4639 1602.8817 1602.8817 1625.4639 55.594925 -101.03316 -41.689385 1724.973 35.280816 26.589804 -40.046683 63.973353 -78.086653 41.559746 41.559746 26.589804 55.594925 + 40 0.12149354 -6.5076059 -4.4000835 2234.612 2345.8821 2337.9552 1513.845 1486.9674 1564.2263 1564.2263 1486.9674 56.775913 -199.0226 -45.302207 1513.845 27.333643 -9.5839331 -16.768077 120.5015 -141.06684 17.802047 17.802047 -9.5839331 56.775913 + 50 0.13053032 -6.5204024 -4.5064456 2097.1985 2270.3252 2180.7348 1458.5923 1396.6336 1507.4402 1507.4402 1396.6336 37.603097 -213.85355 22.82129 1458.5923 47.47204 -20.546302 54.328025 195.36071 -169.33495 -24.624815 -24.624815 -20.546302 37.603097 +Loop time of 0.00293746 on 4 procs for 50 steps with 108 atoms + +Performance: 7353280.980 tau/day, 17021.484 timesteps/s +99.2% 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.00037884 | 0.00040147 | 0.00041828 | 0.0 | 13.67 +Neigh | 0.0010222 | 0.0011142 | 0.0011623 | 0.2 | 37.93 +Comm | 0.00098867 | 0.0010309 | 0.0011341 | 0.2 | 35.09 +Output | 0.00030288 | 0.00031473 | 0.00034979 | 0.0 | 10.71 +Modify | 2.0147e-05 | 2.218e-05 | 2.3691e-05 | 0.0 | 0.76 +Other | | 5.406e-05 | | | 1.84 + +Nlocal: 27 ave 31 max 25 min +Histogram: 2 0 0 1 0 0 0 0 0 1 +Nghost: 429.75 ave 450 max 409 min +Histogram: 1 0 0 0 0 2 0 0 0 1 +Neighs: 727.25 ave 845 max 651 min +Histogram: 1 0 1 1 0 0 0 0 0 1 + +Total # of neighbors = 2909 +Ave neighs/atom = 26.935185 +Neighbor list builds = 50 +Dangerous builds not checked + +Total wall time: 0:00:00 diff --git a/src/EXTRA-COMPUTE/compute_born_matrix.cpp b/src/EXTRA-COMPUTE/compute_born_matrix.cpp new file mode 100644 index 0000000000..2cf0b2decc --- /dev/null +++ b/src/EXTRA-COMPUTE/compute_born_matrix.cpp @@ -0,0 +1,1147 @@ +// clang-format off +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, 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 : Germain Clavier (TUe), Aidan Thompson (Sandia) +------------------------------------------------------------------------- */ + +#include "compute_born_matrix.h" + +#include "angle.h" +#include "atom.h" +#include "atom_vec.h" +#include "bond.h" +#include "comm.h" +#include "compute.h" +#include "dihedral.h" +#include "domain.h" +#include "error.h" +#include "force.h" +#include "improper.h" +#include "kspace.h" +#include "memory.h" +#include "modify.h" +#include "molecule.h" +#include "neigh_list.h" +#include "neigh_request.h" +#include "neighbor.h" +#include "pair.h" +#include "update.h" +#include "universe.h" + +#include +#include + +using namespace LAMMPS_NS; + +#define BIG 1000000000 + +static int constexpr albe[21][2] = { + {0, 0}, // C11 + {1, 1}, // C22 + {2, 2}, // C33 + {1, 2}, // C44 + {0, 2}, // C55 + {0, 1}, // C66 + {0, 1}, // C12 + {0, 2}, // C13 + {0, 3}, // C14 + {0, 4}, // C15 + {0, 5}, // C16 + {1, 2}, // C23 + {1, 3}, // C24 + {1, 4}, // C25 + {1, 5}, // C26 + {2, 3}, // C34 + {2, 4}, // C35 + {2, 5}, // C36 + {3, 4}, // C45 + {3, 5}, // C46 + {4, 5} // C56 +}; + +static int constexpr albemunu[21][4] = { + {0, 0, 0, 0}, // C11 + {1, 1, 1, 1}, // C22 + {2, 2, 2, 2}, // C33 + {1, 2, 1, 2}, // C44 + {0, 2, 0, 2}, // C55 + {0, 1, 0, 1}, // C66 + {0, 0, 1, 1}, // C12 + {0, 0, 2, 2}, // C13 + {0, 0, 1, 2}, // C14 + {0, 0, 0, 2}, // C15 + {0, 0, 0, 1}, // C16 + {1, 1, 2, 2}, // C23 + {1, 1, 1, 2}, // C24 + {1, 1, 0, 2}, // C25 + {1, 1, 0, 1}, // C26 + {2, 2, 1, 2}, // C34 + {2, 2, 0, 2}, // C35 + {2, 2, 0, 1}, // C36 + {1, 2, 0, 2}, // C45 + {1, 2, 0, 1}, // C46 + {0, 1, 0, 2} // C56 +}; + +/* ---------------------------------------------------------------------- */ + +ComputeBornMatrix::ComputeBornMatrix(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg) +{ + if (narg < 3) error->all(FLERR,"Illegal compute born/matrix command"); + + // optional args + + numflag = 0; + numdelta = 0.0; + + int iarg = 3; + while (iarg < narg) { + if (strcmp(arg[iarg],"numdiff") == 0) { + if (iarg+3 > narg) error->all(FLERR,"Illegal compute born/matrix command"); + numflag = 1; + numdelta = utils::numeric(FLERR,arg[iarg+1],false,lmp); + id_virial = utils::strdup(arg[iarg+2]); + int icompute = modify->find_compute(id_virial); + if (icompute < 0) error->all(FLERR,"Could not find compute born/matrix pressure ID"); + compute_virial = modify->compute[icompute]; + if (compute_virial->pressflag == 0) + error->all(FLERR,"Compute born/matrix pressure ID does not compute pressure"); + iarg += 3; + } else error->all(FLERR,"Illegal compute born/matrix command"); + } + + MPI_Comm_rank(world, &me); + + // For now the matrix can be computed as a 21 element vector + + nvalues = 21; + + // Initialize some variables + + values_local = values_global = vector = nullptr; + + // this fix produces a global vector + + memory->create(vector, nvalues, "born_matrix:vector"); + memory->create(values_global, nvalues, "born_matrix:values_global"); + size_vector = nvalues; + + vector_flag = 1; + extvector = 0; + maxatom = 0; + + if (!numflag) { + memory->create(values_local, nvalues, "born_matrix:values_local"); + } else { + + reallocate(); + + // set fixed-point to default = center of cell + + fixedpoint[0] = 0.5 * (domain->boxlo[0] + domain->boxhi[0]); + fixedpoint[1] = 0.5 * (domain->boxlo[1] + domain->boxhi[1]); + fixedpoint[2] = 0.5 * (domain->boxlo[2] + domain->boxhi[2]); + + // define the cartesian indices for each strain (Voigt order) + + dirlist[0][0] = 0; + dirlist[0][1] = 0; + dirlist[1][0] = 1; + dirlist[1][1] = 1; + dirlist[2][0] = 2; + dirlist[2][1] = 2; + + dirlist[3][0] = 1; + dirlist[3][1] = 2; + dirlist[4][0] = 0; + dirlist[4][1] = 2; + dirlist[5][0] = 0; + dirlist[5][1] = 1; + } + +} + +/* ---------------------------------------------------------------------- */ + +ComputeBornMatrix::~ComputeBornMatrix() +{ + memory->destroy(values_global); + memory->destroy(vector); + if (!numflag) { + memory->destroy(values_local); + } else { + memory->destroy(temp_x); + memory->destroy(temp_f); + modify->delete_compute(id_virial); + delete[] id_virial; + } + +} + +/* ---------------------------------------------------------------------- */ + +void ComputeBornMatrix::init() +{ + dt = update->dt; + + pairflag = 0; + bondflag = 0; + angleflag = 0; + dihedflag = 0; + impflag = 0; + kspaceflag = 0; + + if (!numflag) { + if (force->pair) pairflag = 1; + if (force->bond) bondflag = 1; + if (force->angle) angleflag = 1; + if (force->dihedral) dihedflag = 1; + if (force->improper) impflag = 1; + if (force->kspace) kspaceflag = 1; + + // Error check + + if (comm->me == 0) { + if (pairflag && (force->pair->born_matrix_enable == 0)) { + error->all(FLERR, "Pair style does not support compute born/matrix"); + } + // if (bondflag && (force->bond->born_matrix_enable == 0)) { + // error->warning(FLERR, "Bond style does not support compute born/matrix"); + // } + // if (angleflag && (force->angle->born_matrix_enable == 0)) { + // error->warning(FLERR, "Angle style does not support compute born/matrix"); + // } + // if (dihedflag && (force->dihedral->born_matrix_enable == 0)) { + // error->warning(FLERR, "Dihedral style does not support compute born/matrix"); + // } + // if (impflag && (force->improper->born_matrix_enable == 0)) { + // error->warning(FLERR, "Improper style does not support compute born/matrix"); + // } + // if (kspaceflag) { + // error->warning(FLERR, "KSPACE contribution not supported by compute born/matrix"); + // } + } + + // need an occasional half neighbor list + + int irequest = neighbor->request((void *) this); + neighbor->requests[irequest]->pair = 0; + neighbor->requests[irequest]->compute = 1; + neighbor->requests[irequest]->occasional = 1; + + } else { + + // check for virial compute + + int icompute = modify->find_compute(id_virial); + if (icompute < 0) error->all(FLERR, "Virial compute ID for compute born/matrix does not exist"); + compute_virial = modify->compute[icompute]; + + // set up reverse index lookup + + for (int m = 0; m < nvalues; m++) { + int a = albe[m][0]; + int b = albe[m][1]; + revalbe[a][b] = m; + revalbe[b][a] = m; + } + + // voigt3VtoM notation in normal physics sense, + // 3x3 matrix and vector indexing + // i-j: (1-1), (2-2), (3-3), (2-3), (1-3), (1-2) + // voigt3VtoM: 1 2 3 4 5 6 + + voigt3VtoM[0][0]=0; // for 1 + voigt3VtoM[0][1]=0; + voigt3VtoM[1][0]=1; // for 2 + voigt3VtoM[1][1]=1; + voigt3VtoM[2][0]=2; // for 3 + voigt3VtoM[2][1]=2; + voigt3VtoM[3][0]=1; // for 4 + voigt3VtoM[3][1]=2; + voigt3VtoM[4][0]=0; // for 5 + voigt3VtoM[4][1]=2; + voigt3VtoM[5][0]=0; // for 6 + voigt3VtoM[5][1]=1; + + // to convert to vector indexing: + // matrix index to vector index, double -> single index + + voigt3MtoV[0][0]=0; voigt3MtoV[0][1]=5; voigt3MtoV[0][2]=4; + voigt3MtoV[1][0]=5; voigt3MtoV[1][1]=1; voigt3MtoV[1][2]=3; + voigt3MtoV[2][0]=4; voigt3MtoV[2][1]=3; voigt3MtoV[2][2]=2; + + // this is just for the virial. + // since they use the xx,yy,zz,xy,xz,yz + // order not the ordinary voigt + + virialMtoV[0][0]=0; virialMtoV[0][1]=3; virialMtoV[0][2]=4; + virialMtoV[1][0]=3; virialMtoV[1][1]=1; virialMtoV[1][2]=5; + virialMtoV[2][0]=4; virialMtoV[2][1]=5; virialMtoV[2][2]=2; + + // set up 3x3 kronecker deltas + + for(int row = 0; row < NXYZ_VIRIAL; row++) + for(int col = 0; col < NXYZ_VIRIAL; col++) + kronecker[row][col] = 0; + for(int row = 0; row < NXYZ_VIRIAL; row++) + kronecker[row][row] = 1; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputeBornMatrix::init_list(int /* id */, NeighList *ptr) +{ + list = ptr; +} + +/* ---------------------------------------------------------------------- + compute output vector +------------------------------------------------------------------------- */ + +void ComputeBornMatrix::compute_vector() +{ + invoked_array = update->ntimestep; + + if (!numflag) { + + // zero out arrays for one sample + + for (int m = 0; m < nvalues; m++) values_local[m] = 0.0; + + // Compute Born contribution + + if (pairflag) compute_pairs(); + // For now these functions are commented + if (bondflag) compute_bonds(); + if (angleflag) compute_angles(); + if (dihedflag) compute_dihedrals(); + + // sum Born contributions over all procs + + MPI_Allreduce(values_local, values_global, nvalues, MPI_DOUBLE, MPI_SUM, world); + + } else { + + // calculate Born matrix using stress finite differences + + compute_numdiff(); + + } + + for (int m = 0; m < nvalues; m++) { vector[m] = values_global[m]; } + +} + +/* ---------------------------------------------------------------------- + compute Born contribution of local proc +------------------------------------------------------------------------- */ + +void ComputeBornMatrix::compute_pairs() + +{ + int i, j, ii, jj, inum, jnum, itype, jtype; + double rsq, factor_coul, factor_lj; + double dupair, du2pair, rinv; + int *ilist, *jlist, *numneigh, **firstneigh; + + int *type = atom->type; + int *mask = atom->mask; + int nlocal = atom->nlocal; + double *special_coul = force->special_coul; + double *special_lj = force->special_lj; + int newton_pair = force->newton_pair; + + // invoke half neighbor list (will copy or build if necessary) + + neighbor->build_one(list); + + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + + // loop over neighbors of my atoms + + Pair *pair = force->pair; + double **cutsq = force->pair->cutsq; + + // Declares born values + + int a, b, c, d; + double xi1, xi2, xi3; + double fi1, fi2, fi3; + double xj1, xj2, xj3; + double rij[3]; + double pair_pref; + double r2inv; + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + if (!(mask[i] & groupbit)) continue; + + xi1 = atom->x[i][0]; + xi2 = atom->x[i][1]; + xi3 = atom->x[i][2]; + itype = type[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + factor_lj = special_lj[sbmask(j)]; + factor_coul = special_coul[sbmask(j)]; + j &= NEIGHMASK; + + if (!(mask[j] & groupbit)) continue; + + xj1 = atom->x[j][0]; + xj2 = atom->x[j][1]; + xj3 = atom->x[j][2]; + rij[0] = xj1 - xi1; + rij[1] = xj2 - xi2; + rij[2] = xj3 - xi3; + rsq = rij[0] * rij[0] + rij[1] * rij[1] + rij[2] * rij[2]; + r2inv = 1.0 / rsq; + rinv = sqrt(r2inv); + jtype = type[j]; + + if (rsq >= cutsq[itype][jtype]) continue; + + if (newton_pair || j < nlocal) { + + // Add contribution to Born tensor + + pair->born_matrix(i, j, itype, jtype, rsq, factor_coul, factor_lj, dupair, du2pair); + pair_pref = du2pair - dupair * rinv; + + // See albemunu in compute_born_matrix.h for indices order. + + a = 0; + b = 0; + c = 0; + d = 0; + for (int m = 0; m < nvalues; m++) { + a = albemunu[m][0]; + b = albemunu[m][1]; + c = albemunu[m][2]; + d = albemunu[m][3]; + values_local[m] += pair_pref * rij[a] * rij[b] * rij[c] * rij[d] * r2inv; + } + } + } + } +} + +/* ---------------------------------------------------------------------- + compute Born matrix using virial stress finite differences +------------------------------------------------------------------------- */ + +void ComputeBornMatrix::compute_numdiff() +{ + double energy; + int m; + + // grow arrays if necessary + + if (atom->nlocal + atom->nghost > maxatom) reallocate(); + + // store copy of current forces for owned and ghost atoms + + double **x = atom->x; + double **f = atom->f; + int nall = atom->nlocal + atom->nghost; + + for (int i = 0; i < nall; i++) + for (int k = 0; k < 3; k++) { + temp_x[i][k] = x[i][k]; + temp_f[i][k] = f[i][k]; + } + + // loop over 6 strain directions + // compute a finite difference force in each dimension + + int flag, allflag; + + for (int idir = 0; idir < NDIR_VIRIAL; idir++) { + displace_atoms(nall, idir, 1.0); + force_clear(nall); + update_virial(); + for (int jdir = 0; jdir < NDIR_VIRIAL; jdir++) { + m = revalbe[idir][jdir]; + values_global[m] = compute_virial->vector[jdir]; + } + + restore_atoms(nall, idir); + displace_atoms(nall, idir, -1.0); + force_clear(nall); + update_virial(); + for (int jdir = 0; jdir < NDIR_VIRIAL; jdir++) { + m = revalbe[idir][jdir]; + values_global[m] -= compute_virial->vector[jdir]; + } + restore_atoms(nall, idir); + } + + // recompute virial so all virial and energy contributions are as before + // this will possibly break compute stress/atom, need to test + + update_virial(); + + // restore original forces for owned and ghost atoms + + for (int i = 0; i < nall; i++) + for (int k = 0; k < 3; k++) + f[i][k] = temp_f[i][k]; + + double nktv2p = force->nktv2p; + double inv_volume = 1.0 / (domain->xprd * domain->yprd * domain->zprd); + // double denominator = -0.5 / numdelta * inv_volume * nktv2p; + double denominator = -0.5 / numdelta; + + for (int m = 0; m < nvalues; m++) values_global[m] *= denominator; + + virial_addon(); +} + +/* ---------------------------------------------------------------------- + displace position of all owned and ghost atoms +------------------------------------------------------------------------- */ + +void ComputeBornMatrix::displace_atoms(int nall, int idir, double magnitude) +{ + double **x = atom->x; + int k = dirlist[idir][0]; + int l = dirlist[idir][1]; + for (int i = 0; i < nall; i++) + x[i][k] = temp_x[i][k] + numdelta * magnitude * (temp_x[i][l] - fixedpoint[l]); +} + +/* ---------------------------------------------------------------------- + restore position of all owned and ghost atoms +------------------------------------------------------------------------- */ + +void ComputeBornMatrix::restore_atoms(int nall, int idir) +{ + double **x = atom->x; + int k = dirlist[idir][0]; + for (int i = 0; i < nall; i++) { x[i][k] = temp_x[i][k]; } +} + +/* ---------------------------------------------------------------------- + evaluate potential forces and virial + same logic as in Verlet +------------------------------------------------------------------------- */ + +void ComputeBornMatrix::update_virial() +{ + int eflag = 0; + int vflag = 1; + + if (pairflag) force->pair->compute(eflag, vflag); + + if (atom->molecular != Atom::ATOMIC) { + if (force->bond) force->bond->compute(eflag, vflag); + if (force->angle) force->angle->compute(eflag, vflag); + if (force->dihedral) force->dihedral->compute(eflag, vflag); + if (force->improper) force->improper->compute(eflag, vflag); + } + + if (kspaceflag) force->kspace->compute(eflag, vflag); + + compute_virial->compute_vector(); +} + + +/* ---------------------------------------------------------------------- + calculate extra virial terms + following code of Dr. Yubao Zhen + Comp. Phys. Comm. 183 (2012) 261-265 +------------------------------------------------------------------------- */ + +void ComputeBornMatrix::virial_addon() +{ + // compute the contribution due to perturbation + // here the addon parts are put into born + // delta_il sigv_jk + delta_ik sigv_jl + + // delta_jl sigv_ik + delta_jk sigv_il + // Note: in calculation kl is all there from 0 to 6, and ij=(id,jd) + // each time there are six numbers passed for (Dijkl+Djikl) + // and the term I need should be div by 2. + // Job is to arrange the 6 numbers with ij indexing to the 21 element data structure. + // the sigv is the virial stress at current time. It is never touched. + // Note the symmetry of (i-j), (k-n), and (ij, kn) + // so we only need to evaluate 6x6 matrix with symmetry + + int kd, nd, id, jd; + + double* sigv = compute_virial->vector; + + for (int idir = 0; idir < NDIR_VIRIAL; idir++) { + int ijvgt = idir; // this is it. + double addon; + + id = voigt3VtoM[idir][0]; // extract the two indicies composing the voigt reprensentation + jd = voigt3VtoM[idir][1]; + + int SHEAR = 0; + if( id != jd) SHEAR = 1; + + // note BornVec must be cleared before. + + for (int knvgt=ijvgt; knvgt < NDIR_VIRIAL; knvgt++) { + kd= voigt3VtoM[knvgt][0]; + nd= voigt3VtoM[knvgt][1]; + addon = kronecker[id][nd]*sigv[virialMtoV[jd][kd]] + kronecker[id][kd]*sigv[virialMtoV[jd][nd]]; + if(SHEAR) + addon = addon + kronecker[jd][nd]*sigv[virialMtoV[id][kd]] + kronecker[jd][kd]*sigv[virialMtoV[id][nd]]; + } + } +} + +/* ---------------------------------------------------------------------- + clear forces needed +------------------------------------------------------------------------- */ + +void ComputeBornMatrix::force_clear(int nall) +{ + double **forces = atom->f; + size_t nbytes = 3 * sizeof(double) * nall; + if (nbytes) memset(&forces[0][0], 0, nbytes); +} + +/* ---------------------------------------------------------------------- + reallocated local per-atoms arrays +------------------------------------------------------------------------- */ + +void ComputeBornMatrix::reallocate() +{ + memory->destroy(temp_x); + memory->destroy(temp_f); + maxatom = atom->nmax; + memory->create(temp_x, maxatom, 3, "born/matrix:temp_x"); + memory->create(temp_f, maxatom, 3, "born/matrix:temp_f"); +} + +/* ---------------------------------------------------------------------- + memory usage of local atom-based arrays +------------------------------------------------------------------------- */ + +double ComputeBornMatrix::memory_usage() +{ + double bytes = 0.0; + bytes += (double) 2 * maxatom * 3 * sizeof(double); + return bytes; +} + + +/* ---------------------------------------------------------------------- + Count bonds and compute bond info on this proc + only count bond once if newton_bond is off + all atoms in interaction must be in group + all atoms in interaction must be known to proc + if bond is deleted or turned off (type <= 0) + do not count or count contribution + COMMENTED FOR NOW +------------------------------------------------------------------------- */ + +void ComputeBornMatrix::compute_bonds() +{ + /* ---------------------------------------------------------------------- + int i,m,n,nb,atom1,atom2,imol,iatom,btype,ivar; + tagint tagprev; + double dx,dy,dz,rsq; + + double **x = atom->x; + double **v = atom->v; + int *type = atom->type; + tagint *tag = atom->tag; + int *num_bond = atom->num_bond; + tagint **bond_atom = atom->bond_atom; + int **bond_type = atom->bond_type; + int *mask = atom->mask; + + int *molindex = atom->molindex; + int *molatom = atom->molatom; + Molecule **onemols = atom->avec->onemols; + + int nlocal = atom->nlocal; + int newton_bond = force->newton_bond; + int molecular = atom->molecular; + + Bond *bond = force->bond; + + int a,b,c,d; + double rij[3]; + double rinv, r2inv; + double pair_pref, dupair, du2pair; + + // loop over all atoms and their bonds + + m = 0; + while (mnum_bond[iatom]; + } + + for (i = 0; i < nb; i++) { + if (molecular == 1) { + btype = bond_type[atom1][i]; + atom2 = atom->map(bond_atom[atom1][i]); + } else { + tagprev = tag[atom1] - iatom - 1; + btype = onemols[imol]->bond_type[iatom][i]; + atom2 = atom->map(onemols[imol]->bond_atom[iatom][i]+tagprev); + } + + if (atom2 < 0 || !(mask[atom2] & groupbit)) continue; + if (newton_bond == 0 && tag[atom1] > tag[atom2]) continue; + if (btype <= 0) continue; + + dx = x[atom2][0] - x[atom1][0]; + dy = x[atom2][1] - x[atom1][1]; + dz = x[atom2][2] - x[atom1][2]; + domain->minimum_image(dx,dy,dz); + rsq = dx*dx + dy*dy + dz*dz; + rij[0] = dx; + rij[1] = dy; + rij[2] = dz; + r2inv = 1.0/rsq; + rinv = sqrt(r2inv); + + pair_pref = 0.0; + dupair = 0.0; + du2pair = 0.0; + bond->born_matrix(btype,rsq,atom1,atom2,dupair,du2pair); + pair_pref = du2pair - dupair*rinv; + + a = 0; + b = 0; + c = 0; + d = 0; + for (i = 0; i<21; i++) { + a = albemunu[i][0]; + b = albemunu[i][1]; + c = albemunu[i][2]; + d = albemunu[i][3]; + values_local[m+i] += pair_pref*rij[a]*rij[b]*rij[c]*rij[d]*r2inv; + } + } + } + m += 21; + } +------------------------------------------------------------------------- */ +} + +/* ---------------------------------------------------------------------- + count angles and compute angle info on this proc + only count if 2nd atom is the one storing the angle + all atoms in interaction must be in group + all atoms in interaction must be known to proc + if bond is deleted or turned off (type <= 0) + do not count or count contribution + COMMENTED FOR NOW +------------------------------------------------------------------------- */ +void ComputeBornMatrix::compute_angles() +{ + /* ---------------------------------------------------------------------- + int i,m,n,na,atom1,atom2,atom3,imol,iatom,atype,ivar; + tagint tagprev; + double delx1,dely1,delz1,delx2,dely2,delz2; + double rsq1,rsq2,r1,r2,cost; + double r1r2, r1r2inv; + double rsq1inv,rsq2inv,r1inv,r2inv,cinv; + double *ptr; + + double **x = atom->x; + tagint *tag = atom->tag; + int *num_angle = atom->num_angle; + tagint **angle_atom1 = atom->angle_atom1; + tagint **angle_atom2 = atom->angle_atom2; + tagint **angle_atom3 = atom->angle_atom3; + int **angle_type = atom->angle_type; + int *mask = atom->mask; + + int *molindex = atom->molindex; + int *molatom = atom->molatom; + Molecule **onemols = atom->avec->onemols; + + int nlocal = atom->nlocal; + int molecular = atom->molecular; + + // loop over all atoms and their angles + + Angle *angle = force->angle; + + int a,b,c,d,e,f; + double duang, du2ang; + double del1[3], del2[3]; + double dcos[6]; + double d2cos[21]; + double d2lncos[21]; + + // Initializing array for intermediate cos derivatives + // w regard to strain + for (i = 0; i < 6; i++) { + dcos[i] = 0; + } + for (i = 0; i < 21; i++) { + d2cos[i] = 0; + d2lncos[i] = 0; + } + + m = 0; + while (m < nvalues) { + for (atom2 = 0; atom2 < nlocal; atom2++) { + if (!(mask[atom2] & groupbit)) continue; + + if (molecular == 1) na = num_angle[atom2]; + else { + if (molindex[atom2] < 0) continue; + imol = molindex[atom2]; + iatom = molatom[atom2]; + na = onemols[imol]->num_angle[iatom]; + } + + for (i = 0; i < na; i++) { + if (molecular == 1) { + if (tag[atom2] != angle_atom2[atom2][i]) continue; + atype = angle_type[atom2][i]; + atom1 = atom->map(angle_atom1[atom2][i]); + atom3 = atom->map(angle_atom3[atom2][i]); + } else { + if (tag[atom2] != onemols[imol]->angle_atom2[atom2][i]) continue; + atype = onemols[imol]->angle_type[atom2][i]; + tagprev = tag[atom2] - iatom - 1; + atom1 = atom->map(onemols[imol]->angle_atom1[atom2][i]+tagprev); + atom3 = atom->map(onemols[imol]->angle_atom3[atom2][i]+tagprev); + } + + if (atom1 < 0 || !(mask[atom1] & groupbit)) continue; + if (atom3 < 0 || !(mask[atom3] & groupbit)) continue; + if (atype <= 0) continue; + + delx1 = x[atom1][0] - x[atom2][0]; + dely1 = x[atom1][1] - x[atom2][1]; + delz1 = x[atom1][2] - x[atom2][2]; + domain->minimum_image(delx1,dely1,delz1); + del1[0] = delx1; + del1[1] = dely1; + del1[2] = delz1; + + rsq1 = delx1*delx1 + dely1*dely1 + delz1*delz1; + rsq1inv = 1.0/rsq1; + r1 = sqrt(rsq1); + r1inv = 1.0/r1; + + delx2 = x[atom3][0] - x[atom2][0]; + dely2 = x[atom3][1] - x[atom2][1]; + delz2 = x[atom3][2] - x[atom2][2]; + domain->minimum_image(delx2,dely2,delz2); + del2[0] = delx2; + del2[1] = dely2; + del2[2] = delz2; + + rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2; + rsq2inv = 1.0/rsq2; + r2 = sqrt(rsq2); + r2inv = 1.0/r2; + + r1r2 = delx1*delx2 + dely1*dely2 + delz1*delz2; + r1r2inv = 1/r1r2; + // cost = cosine of angle + + cost = delx1*delx2 + dely1*dely2 + delz1*delz2; + cost /= r1*r2; + if (cost > 1.0) cost = 1.0; + if (cost < -1.0) cost = -1.0; + cinv = 1.0/cost; + + // The method must return derivative with regards + // to cos(theta)! + // Use the chain rule if needed: + // dU(t)/de = dt/dcos(t)*dU(t)/dt*dcos(t)/de + // with dt/dcos(t) = -1/sin(t) + angle->born_matrix(atype,atom1,atom2,atom3,duang,du2ang); + + // Voigt notation + // 1 = 11, 2 = 22, 3 = 33 + // 4 = 23, 5 = 13, 6 = 12 + a = 0; + b = 0; + c = 0; + d = 0; + for (i = 0; i<6; i++) { + a = albe[i][0]; + b = albe[i][1]; + dcos[i] = cost*(del1[a]*del2[b]+del1[b]*del2[a]*r1r2inv - + del1[a]*del1[b]*rsq1inv - del2[a]*del2[b]*rsq2inv); + } + for (i = 0; i<21; i++) { + a = albemunu[i][0]; + b = albemunu[i][1]; + c = albemunu[i][2]; + d = albemunu[i][3]; + e = albe[i][0]; + f = albe[i][1]; + d2lncos[i] = 2*(del1[a]*del1[b]*del1[c]*del1[d]*rsq1inv*rsq1inv + + del2[a]*del2[b]*del2[c]*del2[d]*rsq2inv*rsq2inv) - + (del1[a]*del2[b]+del1[b]*del2[a]) * + (del1[c]*del2[d]+del1[d]*del2[c]) * + r1r2inv*r1r2inv; + d2cos[i] = cost*d2lncos[i] + dcos[e]*dcos[f]*cinv; + values_local[m+i] += duang*d2cos[i] + du2ang*dcos[e]*dcos[f]; + } + } + } + m+=21; + } +------------------------------------------------------------------------- */ +} + +/* ---------------------------------------------------------------------- + count dihedrals on this proc + only count if 2nd atom is the one storing the dihedral + all atoms in interaction must be in group + all atoms in interaction must be known to proc + if flag is set, compute requested info about dihedral + COMMENTED FOR NOW +------------------------------------------------------------------------- */ + +void ComputeBornMatrix::compute_dihedrals() +{ + /* ---------------------------------------------------------------------- + int i,m,n,nd,atom1,atom2,atom3,atom4,imol,iatom,dtype,ivar; + tagint tagprev; + double vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z,vb2xm,vb2ym,vb2zm; + double ax,ay,az,bx,by,bz,rasq,rbsq,rgsq,rg,ra2inv,rb2inv,rabinv; + double si,co,phi; + double *ptr; + + double **x = atom->x; + tagint *tag = atom->tag; + int *num_dihedral = atom->num_dihedral; + tagint **dihedral_atom1 = atom->dihedral_atom1; + tagint **dihedral_atom2 = atom->dihedral_atom2; + tagint **dihedral_atom3 = atom->dihedral_atom3; + tagint **dihedral_atom4 = atom->dihedral_atom4; + int *mask = atom->mask; + + int *molindex = atom->molindex; + int *molatom = atom->molatom; + Molecule **onemols = atom->avec->onemols; + + int nlocal = atom->nlocal; + int molecular = atom->molecular; + + // loop over all atoms and their dihedrals + + Dihedral *dihedral = force->dihedral; + + double dudih,du2dih; + int a,b,c,d,e,f; + double b1sq; + double b2sq; + double b3sq; + double b1b2; + double b1b3; + double b2b3; + double b1[3]; + double b2[3]; + double b3[3]; + + // Actually derivatives of the square of the + // vectors dot product. + double dmn[6]; + double dmm[6]; + double dnn[6]; + double d2mn[21]; + double d2mm[21]; + double d2nn[21]; + + double dcos[6]; + double d2cos[21]; + + for (i = 0; i < 6; i++) { + dmn[i] =0; + dmm[i] = 0; + dnn[i] = 0; + dcos[i] = 0; + } + for (i = 0; i < 21; i++) { + d2mn[i] = 0; + d2mm[i] = 0; + d2nn[i] = 0; + d2cos[i] = 0; + } + + m = 0; + while (m < nvalues) { + for (atom2 = 0; atom2 < nlocal; atom2++) { + if (!(mask[atom2] & groupbit)) continue; + + if (molecular == 1) nd = num_dihedral[atom2]; + else { + if (molindex[atom2] < 0) continue; + imol = molindex[atom2]; + iatom = molatom[atom2]; + nd = onemols[imol]->num_dihedral[iatom]; + } + + for (i = 0; i < nd; i++) { + if (molecular == 1) { + if (tag[atom2] != dihedral_atom2[atom2][i]) continue; + atom1 = atom->map(dihedral_atom1[atom2][i]); + atom3 = atom->map(dihedral_atom3[atom2][i]); + atom4 = atom->map(dihedral_atom4[atom2][i]); + } else { + if (tag[atom2] != onemols[imol]->dihedral_atom2[atom2][i]) continue; + tagprev = tag[atom2] - iatom - 1; + atom1 = atom->map(onemols[imol]->dihedral_atom1[atom2][i]+tagprev); + atom3 = atom->map(onemols[imol]->dihedral_atom3[atom2][i]+tagprev); + atom4 = atom->map(onemols[imol]->dihedral_atom4[atom2][i]+tagprev); + } + + if (atom1 < 0 || !(mask[atom1] & groupbit)) continue; + if (atom3 < 0 || !(mask[atom3] & groupbit)) continue; + if (atom4 < 0 || !(mask[atom4] & groupbit)) continue; + + // phi calculation from dihedral style harmonic + + // The method must return derivative with regards + // to cos(phi)! + // Use the chain rule if needed: + // dU(t)/de = dt/dcos(t)*dU(t)/dt*dcos(t)/de + // with dt/dcos(t) = -1/sin(t) + + dihedral->born_matrix(nd,atom1,atom2,atom3,atom4,dudih,du2dih); + + vb1x = x[atom1][0] - x[atom2][0]; + vb1y = x[atom1][1] - x[atom2][1]; + vb1z = x[atom1][2] - x[atom2][2]; + domain->minimum_image(vb1x,vb1y,vb1z); + b1[0] = vb1x; + b1[1] = vb1y; + b1[2] = vb1z; + b1sq = b1[0]*b1[0]+b1[1]*b1[1]+b1[2]*b1[2]; + + vb2x = x[atom3][0] - x[atom2][0]; + vb2y = x[atom3][1] - x[atom2][1]; + vb2z = x[atom3][2] - x[atom2][2]; + domain->minimum_image(vb2x,vb2y,vb2z); + b2[0] = vb2x; + b2[1] = vb2y; + b2[2] = vb2z; + b2sq = b2[0]*b2[0]+b2[1]*b2[1]+b2[2]*b2[2]; + + vb2xm = -vb2x; + vb2ym = -vb2y; + vb2zm = -vb2z; + domain->minimum_image(vb2xm,vb2ym,vb2zm); + + vb3x = x[atom4][0] - x[atom3][0]; + vb3y = x[atom4][1] - x[atom3][1]; + vb3z = x[atom4][2] - x[atom3][2]; + domain->minimum_image(vb3x,vb3y,vb3z); + b3[0] = vb3x; + b3[1] = vb3y; + b3[2] = vb3z; + b3sq = b3[0]*b3[0]+b3[1]*b3[1]+b3[2]*b3[2]; + + b1b2 = b1[0]*b2[0]+b1[1]*b2[1]+b1[2]*b2[2]; + b1b3 = b1[0]*b3[0]+b1[1]*b3[1]+b1[2]*b3[2]; + b2b3 = b2[0]*b3[0]+b2[1]*b3[1]+b2[2]*b3[2]; + + ax = vb1y*vb2zm - vb1z*vb2ym; + ay = vb1z*vb2xm - vb1x*vb2zm; + az = vb1x*vb2ym - vb1y*vb2xm; + bx = vb3y*vb2zm - vb3z*vb2ym; + by = vb3z*vb2xm - vb3x*vb2zm; + bz = vb3x*vb2ym - vb3y*vb2xm; + + rasq = ax*ax + ay*ay + az*az; + rbsq = bx*bx + by*by + bz*bz; + rgsq = vb2xm*vb2xm + vb2ym*vb2ym + vb2zm*vb2zm; + rg = sqrt(rgsq); + + ra2inv = rb2inv = 0.0; + if (rasq > 0) ra2inv = 1.0/rasq; + if (rbsq > 0) rb2inv = 1.0/rbsq; + rabinv = sqrt(ra2inv*rb2inv); + + co = (ax*bx + ay*by + az*bz)*rabinv; + si = rg*rabinv*(ax*vb3x + ay*vb3y + az*vb3z); + + if (co > 1.0) co = 1.0; + if (co < -1.0) co = -1.0; + phi = atan2(si,co); + + // above a and b are m and n vectors + // here they are integers indices + a = 0; + b = 0; + c = 0; + d = 0; + e = 0; + f = 0; + for (i = 0; i<6; i++) { + a = albe[i][0]; + b = albe[i][1]; + dmm[i] = 2*(b2sq*b1[a]*b1[b]+b1sq*b2[a]*b2[b] - b1b2*(b1[a]*b2[b]+b1[b]*b2[a])); + dnn[i] = 2*(b3sq*b2[a]*b2[b]+b2sq*b3[a]*b3[b] - b2b3*(b2[a]*b3[b]+b2[b]*b3[a])); + dmn[i] = b1b2*(b2[a]*b3[b]+b2[b]*b3[a]) + b2b3*(b1[a]*b2[b]+b1[b]*b2[a]) + - 2*(b1b3*b2[a]*b2[b]) - b2sq*(b1[a]*b3[b]+b1[b]*b3[a]); + dcos[i] = co*(rabinv*rabinv*dmn[i] - ra2inv*dmm[i] - rb2inv*dnn[i])/2.; + } + for (i = 0; i<21; i++) { + a = albemunu[i][0]; + b = albemunu[i][1]; + c = albemunu[i][2]; + d = albemunu[i][3]; + e = albe[i][0]; + f = albe[i][1]; + d2mm[i] = 4*(b1[a]*b1[b]*b2[c]*b2[d] + b1[c]*b1[d]*b2[a]*b2[b]) + - 8*(b1[a]*b2[b]+b1[b]*b2[a])*(b1[c]*b2[d]+b1[d]*b2[c]); + d2nn[i] = 4*(b2[a]*b2[b]*b3[c]*b3[d] + b2[c]*b2[d]*b3[a]*b3[b]) + - 8*(b2[a]*b3[b]+b2[b]*b3[a])*(b2[c]*b3[d]+b2[d]*b3[c]); + d2mn[i] = (b1[a]*b2[b]+b1[b]*b2[a])*(b2[c]*b3[d]+b2[d]*b3[c]) + + (b2[a]*b3[b]+b2[b]*b3[a])*(b1[c]*b2[d]+b1[d]*b2[d]) + - (b1[a]*b3[b]+b1[b]*b3[a])*(b2[c]*b2[d]+b2[c]*b2[d]) + - (b1[c]*b3[d]+b1[d]*b3[c])*(b2[a]*b2[b]+b2[a]*b2[b]); + d2cos[i] = co/2.*( + rabinv*rabinv*d2mn[i] + - rabinv*rabinv*rabinv*rabinv*dmn[e]*dmn[f] + + ra2inv*ra2inv*dmm[e]*dmm[f] + - ra2inv*d2mm[i] + + rb2inv*rb2inv*dnn[e]*dnn[f] + - rb2inv*d2nn[i] ); + values_local[m+i] += dudih*d2cos[i] + du2dih*dcos[e]*dcos[f]; + } + } + } + m+=21; + } +------------------------------------------------------------------------- */ +} + diff --git a/src/EXTRA-COMPUTE/compute_born_matrix.h b/src/EXTRA-COMPUTE/compute_born_matrix.h new file mode 100644 index 0000000000..cc54048ce2 --- /dev/null +++ b/src/EXTRA-COMPUTE/compute_born_matrix.h @@ -0,0 +1,103 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://www.lammps.org/, 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 : Germain Clavier (TUe) + --------------------------------------------------------------------------*/ + +#ifdef COMPUTE_CLASS +// clang-format off +ComputeStyle(born/matrix,ComputeBornMatrix); +// clang-format on +#else + +#ifndef LMP_COMPUTE_BORN_MATRIX_H +#define LMP_COMPUTE_BORN_MATRIX_H + +#include "compute.h" + +namespace LAMMPS_NS { + +class ComputeBornMatrix : public Compute { + public: + ComputeBornMatrix(class LAMMPS *, int, char **); + virtual ~ComputeBornMatrix() override; + void init() override; + void init_list(int, class NeighList *) override; + void compute_vector() override; + double memory_usage() override; + + private: + + // Born matrix contributions + + void compute_pairs(); // pair and manybody + void compute_bonds(); // bonds + void compute_angles(); // angles + void compute_dihedrals(); // dihedrals + void compute_numdiff(); // stress virial finite differences + void displace_atoms(int, int, double); // displace atoms + void force_clear(int); // zero out force array + void update_virial(); // recalculate the virial + void restore_atoms(int, int); // restore atom positions + void virial_addon(); // restore atom positions + void reallocate(); // grow the atom arrays + + int me; // process rank + int nvalues; // length of elastic tensor + int numflag; // 1 if using finite differences + double numdelta; // size of finite strain + int maxatom; // allocated size of atom arrays + + int pairflag, bondflag, angleflag; + int dihedflag, impflag, kspaceflag; + + double *values_local,*values_global; + double pos,pos1,dt,nktv2p,ftm2v; + class NeighList *list; + + char *id_virial; // name of virial compute + class Compute *compute_virial; // pointer to virial compute + + static constexpr int NDIR_VIRIAL = 6; // dimension of virial and strain vectors + static constexpr int NXYZ_VIRIAL = 3; // number of Cartesian coordinates + int revalbe[NDIR_VIRIAL][NDIR_VIRIAL]; + int voigt3VtoM[NDIR_VIRIAL][2]; + int voigt3MtoV[NXYZ_VIRIAL][NXYZ_VIRIAL]; + int virialMtoV[NXYZ_VIRIAL][NXYZ_VIRIAL]; + int kronecker[NXYZ_VIRIAL][NXYZ_VIRIAL]; + double **temp_x; // original coords + double **temp_f; // original forces + double fixedpoint[NXYZ_VIRIAL]; // displacement field origin + int dirlist[NDIR_VIRIAL][2]; // strain cartesian indices +}; + +} + +#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: ... style does not support compute born/matrix + +Some component of the force field (pair, bond, angle...) does not provide +a function to return the Born term contribution. +*/ + diff --git a/src/pair.cpp b/src/pair.cpp index f88c4e0972..23e6ad5bf2 100644 --- a/src/pair.cpp +++ b/src/pair.cpp @@ -58,6 +58,7 @@ Pair::Pair(LAMMPS *lmp) : Pointers(lmp) comm_forward = comm_reverse = comm_reverse_off = 0; single_enable = 1; + born_matrix_enable = 0; single_hessian_enable = 0; restartinfo = 1; respa_enable = 0; diff --git a/src/pair.h b/src/pair.h index a61af0c1da..addf4a2627 100644 --- a/src/pair.h +++ b/src/pair.h @@ -52,6 +52,7 @@ class Pair : protected Pointers { int single_enable; // 1 if single() routine exists int single_hessian_enable; // 1 if single_hessian() routine exists + int born_matrix_enable; // 1 if born_matrix() routine exists int restartinfo; // 1 if pair style writes restart info int respa_enable; // 1 if inner/middle/outer rRESPA routines int one_coeff; // 1 if allows only one coeff * * call @@ -168,6 +169,13 @@ class Pair : protected Pointers { return 0.0; } + virtual void born_matrix(int /*i*/, int /*j*/, int /*itype*/, int /*jtype*/, double /*rsq*/, + double /*factor_coul*/, double /*factor_lj*/, double& du, double& du2) + { + du = 0.0; + du2 = 0.0; + } + virtual void settings(int, char **) = 0; virtual void coeff(int, char **) = 0; diff --git a/src/pair_lj_cut.cpp b/src/pair_lj_cut.cpp index a9d45b9007..983798d3dd 100644 --- a/src/pair_lj_cut.cpp +++ b/src/pair_lj_cut.cpp @@ -40,6 +40,7 @@ using namespace MathConst; PairLJCut::PairLJCut(LAMMPS *lmp) : Pair(lmp) { respa_enable = 1; + born_matrix_enable = 1; writedata = 1; } @@ -680,6 +681,28 @@ double PairLJCut::single(int /*i*/, int /*j*/, int itype, int jtype, double rsq, /* ---------------------------------------------------------------------- */ +void PairLJCut::born_matrix(int /*i*/, int /*j*/, int itype, int jtype, double rsq, + double /*factor_coul*/, double factor_lj, + double &dupair, double &du2pair) +{ + double rinv,r2inv,r6inv,du,du2; + + r2inv = 1.0/rsq; + rinv = sqrt(r2inv); + r6inv = r2inv*r2inv*r2inv; + + // Reminder: lj1 = 48*e*s^12, lj2 = 24*e*s^6 + // so dupair = -forcelj/r = -fforce*r (forcelj from single method) + + du = r6inv * rinv * (lj2[itype][jtype] - lj1[itype][jtype]*r6inv); + du2 = r6inv * r2inv * (13*lj1[itype][jtype]*r6inv - 7*lj2[itype][jtype]); + + dupair = factor_lj*du; + du2pair = factor_lj*du2; +} + +/* ---------------------------------------------------------------------- */ + void *PairLJCut::extract(const char *str, int &dim) { dim = 2; diff --git a/src/pair_lj_cut.h b/src/pair_lj_cut.h index 8ca9a14620..76bd29eb05 100644 --- a/src/pair_lj_cut.h +++ b/src/pair_lj_cut.h @@ -40,6 +40,7 @@ class PairLJCut : public Pair { void write_data(FILE *) override; void write_data_all(FILE *) override; double single(int, int, int, int, double, double, double, double &) override; + void born_matrix(int, int, int, int, double, double, double, double &, double &) override; void *extract(const char *, int &) override; void compute_inner() override;