diff --git a/doc/src/Commands_compute.rst b/doc/src/Commands_compute.rst index 9b04bbff92..61c5e83eda 100644 --- a/doc/src/Commands_compute.rst +++ b/doc/src/Commands_compute.rst @@ -33,6 +33,7 @@ KOKKOS, o = OPENMP, t = OPT. * :doc:`body/local ` * :doc:`bond ` * :doc:`bond/local ` + * :doc:`born/matrix ` * :doc:`centro/atom ` * :doc:`centroid/stress/atom ` * :doc:`chunk/atom ` diff --git a/doc/src/Developer_utils.rst b/doc/src/Developer_utils.rst index aceecd7a46..d7b1f077dd 100644 --- a/doc/src/Developer_utils.rst +++ b/doc/src/Developer_utils.rst @@ -214,6 +214,9 @@ Convenience functions .. doxygenfunction:: errorurl :project: progguide +.. doxygenfunction:: missing_cmd_args + :project: progguide + .. doxygenfunction:: flush_buffers(LAMMPS *lmp) :project: progguide diff --git a/doc/src/Howto_elastic.rst b/doc/src/Howto_elastic.rst index 4870942458..d77ff3e778 100644 --- a/doc/src/Howto_elastic.rst +++ b/doc/src/Howto_elastic.rst @@ -18,23 +18,52 @@ At zero temperature, it is easy to estimate these derivatives by deforming the simulation box in one of the six directions using the :doc:`change_box ` command and measuring the change in the stress tensor. A general-purpose script that does this is given in the -examples/elastic directory described on the :doc:`Examples ` +examples/ELASTIC directory described on the :doc:`Examples ` doc page. Calculating elastic constants at finite temperature is more challenging, because it is necessary to run a simulation that performs -time averages of differential properties. One way to do this is to -measure the change in average stress tensor in an NVT simulations when +time averages of differential properties. There are at least +3 ways to do this in LAMMPS. The most reliable way to do this is +by exploiting the relationship between elastic constants, stress +fluctuations, and the Born matrix, the second derivatives of energy +w.r.t. strain :ref:`(Ray) `. +The Born matrix calculation has been enabled by +the :doc:`compute born/matrix ` command, +which works for any bonded or non-bonded potential in LAMMPS. +The most expensive part of the calculation is the sampling of +the stress fluctuations. Several examples of this method are +provided in the examples/ELASTIC_T/BORN_MATRIX directory +described on the :doc:`Examples ` doc page. + +A second way is to measure +the change in average stress tensor in an NVT simulations when the cell volume undergoes a finite deformation. In order to balance the systematic and statistical errors in this method, the magnitude of the deformation must be chosen judiciously, and care must be taken to fully equilibrate the deformed cell before sampling the stress -tensor. Another approach is to sample the triclinic cell fluctuations +tensor. An example of this method is provided in the +examples/ELASTIC_T/DEFORMATION directory +described on the :doc:`Examples ` doc page. + +Another approach is to sample the triclinic cell fluctuations that occur in an NPT simulation. This method can also be slow to -converge and requires careful post-processing :ref:`(Shinoda) ` +converge and requires careful post-processing :ref:`(Shinoda) `. +We do not provide an example of this method. + +A nice review of the advantages and disadvantages of all of these methods +is provided in the paper by Clavier et al. :ref:`(Clavier) `. ---------- +.. _Ray: + +**(Ray)** J. R. Ray and A. Rahman, J Chem Phys, 80, 4423 (1984). + .. _Shinoda1: **(Shinoda)** Shinoda, Shiga, and Mikami, Phys Rev B, 69, 134103 (2004). + +.. _Clavier: + +**(Clavier)** G. Clavier, N. Desbiens, E. Bourasseau, V. Lachet, N. Brusselle-Dupend and B. Rousseau, Mol Sim, 43, 1413 (2017). diff --git a/doc/src/compute.rst b/doc/src/compute.rst index 6da8455165..6fdedbbb95 100644 --- a/doc/src/compute.rst +++ b/doc/src/compute.rst @@ -179,6 +179,7 @@ The individual style names on the :doc:`Commands compute ` pag * :doc:`body/local ` - attributes of body sub-particles * :doc:`bond ` - energy of each bond sub-style * :doc:`bond/local ` - distance and energy of each bond +* :doc:`born/matrix ` - second derivative or potential with respect to strain * :doc:`centro/atom ` - centro-symmetry parameter for each atom * :doc:`centroid/stress/atom ` - centroid based stress tensor for each atom * :doc:`chunk/atom ` - assign chunk IDs to each atom diff --git a/doc/src/compute_born_matrix.rst b/doc/src/compute_born_matrix.rst new file mode 100644 index 0000000000..84d002a621 --- /dev/null +++ b/doc/src/compute_born_matrix.rst @@ -0,0 +1,213 @@ +.. index:: compute born/matrix + +compute born/matrix command +=========================== + +Syntax +"""""" + +.. parsed-literal:: + + compute ID group-ID born/matrix keyword value ... + +* ID, group-ID are documented in :doc:`compute ` command +* born/matrix = style name of this compute command +* zero or more keyword/value pairs may be appended + + .. parsed-literal:: + + keyword = *numdiff* + *numdiff* values = delta virial-ID + delta = magnitude of strain (dimensionless) + virial-ID = ID of pressure compute for virial (string) + +Examples +"""""""" + +.. code-block:: LAMMPS + + compute 1 all born/matrix + compute 1 all born/matrix bond angle + compute 1 all born/matrix numdiff 1.0e-4 myvirial + +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` w.r.t. 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, as is the matrix of second derivatives +of potential energy w.r.t strain, +whose 21 independent elements are output in this order: + +.. 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 element +:math:`k` in the global vector output by this compute. Each term comes from the sum +of the derivatives of every contribution to the potential energy +in the system as explained in :ref:`(VanWorkum) +`. + +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. Several examples of this method are +provided in the examples/ELASTIC_T/BORN_MATRIX directory +described on the :doc:`Examples ` doc page. + +NOTE: In the above :math:`C_{i,j}` computation, the fluctuation +term involving the virial stress tensor :math:`\sigma` is the +covariance between each elements. In a +solid the stress fluctuations can vary rapidly, while average +fluctuations can be slow to converge. +A detailed analysis of the convergence rate of all the terms in +the elastic tensor +is provided in the paper by Clavier et al. :ref:`(Clavier) `. + +Two different computation methods for the Born matrix are implemented in this +compute and are mutually exclusive. + +The first one is a direct computation from the analytical formula from the +different terms of the potential used for the simulations :ref:`(VanWorkum) +`. However, the implementation of such derivations must be done +for every potential form. This has not been done yet and can be very +complicated for complex potentials. At the moment a warning message is +displayed for every term that is not supporting the compute at the moment. +This method is the default for now. + +The second method uses finite differences of energy to numerically approximate +the second derivatives :ref:`(Zhen) `. This is useful when using +interaction styles for which the analytical second derivatives have not been +implemented. In this cases, the compute applies linear strain fields of +magnitude *delta* to all the atoms relative to a point at the center of the +box. The strain fields are in six different directions, corresponding to the +six Cartesian components of the stress tensor defined by LAMMPS. For each +direction it applies the strain field in both the positive and negative senses, +and the new stress virial tensor of the entire system is calculated after each. +The difference in these two virials divided by two times *delta*, approximates +the corresponding components of the second derivative, after applying a +suitable unit conversion. + +.. note:: + + It is important to choose a suitable value for delta, the magnitude of + strains that are used to generate finite difference + approximations to the exact virial stress. For typical systems, a value in + the range of 1 part in 1e5 to 1e6 will be sufficient. + However, the best value will depend on a multitude of factors + including the stiffness of the interatomic potential, the thermodynamic + state of the material being probed, and so on. The only way to be sure + that you have made a good choice is to do a sensitivity study on a + representative atomic configuration, sweeping over a wide range of + values of delta. If delta is too small, the output values will vary + erratically due to truncation effects. If delta is increased beyond a + certain point, the output values will start to vary smoothly with + delta, due to growing contributions from higher order derivatives. In + between these two limits, the numerical virial values should be largely + independent of delta. + +The keyword requires the additional arguments *delta* and *virial-ID*. +*delta* gives the size of the applied strains. *virial-ID* gives +the ID string of the pressure compute that provides the virial stress tensor, +requiring that it use the virial keyword e.g. + +.. code-block:: LAMMPS + + compute myvirial all pressure NULL virial + compute 1 all born/matrix numdiff 1.0e-4 myvirial + +**Output info:** + +This compute calculates a global vector with 21 values that are +the second derivatives of the potential energy w.r.t. strain. +The values are in energy units. +The values are ordered as explained above. These values can be used +by any command that uses 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. 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 LAMMPS interaction styles implement the *born_matrix* method +giving first and second order derivatives and LAMMPS will exit with an error if +this compute is used with such interactions unless the *numdiff* option is +also used. The *numdiff* option cannot be used with any other keyword. In this +situation, LAMMPS will also exit with an error. + +Default +""""""" + +none + +---------- + +.. _VanWorkum: + +**(Van Workum)** K. Van Workum et al., J. Chem. Phys. 125 144506 (2006) + +.. _Clavier2: + +**(Clavier)** G. Clavier, N. Desbiens, E. Bourasseau, V. Lachet, N. Brusselle-Dupend and B. Rousseau, Mol Sim, 43, 1413 (2017). + +.. _Zhen: + +**(Zhen)** Y. Zhen, C. Chu, Computer Physics Communications 183(2012)261-265 diff --git a/doc/src/compute_temp_profile.rst b/doc/src/compute_temp_profile.rst index b70206781f..e5830f5cf7 100644 --- a/doc/src/compute_temp_profile.rst +++ b/doc/src/compute_temp_profile.rst @@ -76,21 +76,28 @@ velocity for each atom. Note that if there is only one atom in the bin, its thermal velocity will thus be 0.0. After the spatially-averaged velocity field has been subtracted from -each atom, the temperature is calculated by the formula KE = (dim\*N -- dim\*Nx\*Ny\*Nz) k T/2, where KE = total kinetic energy of the group of -atoms (sum of 1/2 m v\^2), dim = 2 or 3 = dimensionality of the -simulation, N = number of atoms in the group, k = Boltzmann constant, -and T = temperature. The dim\*Nx\*Ny\*Nz term are degrees of freedom -subtracted to adjust for the removal of the center-of-mass velocity in -each of Nx\*Ny\*Nz bins, as discussed in the :ref:`(Evans) ` paper. +each atom, the temperature is calculated by the formula +*KE* = (*dim\*N* - *Ns\*Nx\*Ny\*Nz* - *extra* ) *k* *T*/2, where *KE* = total +kinetic energy of the group of atoms (sum of 1/2 *m* *v*\^2), *dim* = 2 +or 3 = dimensionality of the simulation, *Ns* = 0, 1, 2 or 3 for +streaming velocity subtracted in 0, 1, 2 or 3 dimensions, *extra* = extra +degrees-of-freedom, *N* = number of atoms in the group, *k* = Boltzmann +constant, and *T* = temperature. The *Ns\*Nx\*Ny\*Nz* term is degrees +of freedom subtracted to adjust for the removal of the center-of-mass +velocity in each direction of the *Nx\*Ny\*Nz* bins, as discussed in the +:ref:`(Evans) ` paper. The extra term defaults to (*dim* - *Ns*) +and accounts for overall conservation of center-of-mass velocity across +the group in directions where streaming velocity is *not* subtracted. This +can be altered using the *extra* option of the +:doc:`compute_modify ` command. If the *out* keyword is used with a *tensor* value, which is the default, a kinetic energy tensor, stored as a 6-element vector, is also calculated by this compute for use in the computation of a pressure tensor. The formula for the components of the tensor is the -same as the above formula, except that v\^2 is replaced by vx\*vy for -the xy component, etc. The 6 components of the vector are ordered xx, -yy, zz, xy, xz, yz. +same as the above formula, except that *v*\^2 is replaced by *vx\*vy* for +the xy component, etc. The 6 components of the vector are ordered *xx, +yy, zz, xy, xz, yz.* If the *out* keyword is used with a *bin* value, the count of atoms and computed temperature for each bin are stored for output, as an @@ -123,10 +130,20 @@ needed, the subtracted degrees-of-freedom can be altered using the .. note:: When using the *out* keyword with a value of *bin*, the - calculated temperature for each bin does not include the - degrees-of-freedom adjustment described in the preceding paragraph, - for fixes that constrain molecular motion. It does include the - adjustment due to the *extra* option, which is applied to each bin. + calculated temperature for each bin includes the degrees-of-freedom + adjustment described in the preceding paragraph for fixes that + constrain molecular motion, as well as the adjustment due to + the *extra* option (which defaults to *dim* - *Ns* as described above), + by fractionally applying them based on the fraction of atoms in each + bin. As a result, the bin degrees-of-freedom summed over all bins exactly + equals the degrees-of-freedom used in the scalar temperature calculation, + :math:`\Sigma N_{DOF_i} = N_{DOF}` and the corresponding relation for temperature + is also satisfied :math:`\Sigma N_{DOF_i} T_i = N_{DOF} T`. + These relations will breakdown in cases where the adjustment + exceeds the actual number of degrees-of-freedom in a bin. This could happen + if a bin is empty or in situations where rigid molecules + are non-uniformly distributed, in which case the reported + temperature within a bin may not be accurate. See the :doc:`Howto thermostat ` page for a discussion of different ways to compute temperature and perform diff --git a/doc/utils/sphinx-config/false_positives.txt b/doc/utils/sphinx-config/false_positives.txt index ca62bd9c41..ee183a5877 100644 --- a/doc/utils/sphinx-config/false_positives.txt +++ b/doc/utils/sphinx-config/false_positives.txt @@ -3599,6 +3599,7 @@ Voronoi VORONOI Vorselaars Voth +Voyiatzis vpz vratio Vries @@ -3664,7 +3665,7 @@ Wittmaack wn Wolde workflow -workflows +Workum Worley Wriggers Wuppertal diff --git a/examples/ELASTIC_T/BORN_MATRIX/Argon/Analytical/born_matrix.out b/examples/ELASTIC_T/BORN_MATRIX/Argon/Analytical/born_matrix.out new file mode 100644 index 0000000000..e6854d7bfe --- /dev/null +++ b/examples/ELASTIC_T/BORN_MATRIX/Argon/Analytical/born_matrix.out @@ -0,0 +1,7 @@ +# Cij Matrix from post process computation + 3.36316 1.87373 1.87607 -0.00346 -0.00172 -0.00104 + 1.87373 3.36170 1.87425 0.00443 0.00033 0.00014 + 1.87607 1.87425 3.36573 0.00143 0.00155 0.00127 +-0.00346 0.00443 0.00143 1.87425 0.00127 0.00033 +-0.00172 0.00033 0.00155 0.00127 1.87607 -0.00346 +-0.00104 0.00014 0.00127 0.00033 -0.00346 1.87373 diff --git a/examples/ELASTIC_T/BORN_MATRIX/Argon/Analytical/compute_born.py b/examples/ELASTIC_T/BORN_MATRIX/Argon/Analytical/compute_born.py new file mode 100644 index 0000000000..0f3265faeb --- /dev/null +++ b/examples/ELASTIC_T/BORN_MATRIX/Argon/Analytical/compute_born.py @@ -0,0 +1,118 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +import sys +import numpy as np + +def reduce_Born(Cf): + C = np.zeros((6,6), dtype=np.float64) + C[0,0] = Cf[0] + C[1,1] = Cf[1] + C[2,2] = Cf[2] + C[3,3] = Cf[3] + C[4,4] = Cf[4] + C[5,5] = Cf[5] + C[0,1] = Cf[6] + C[0,2] = Cf[7] + C[0,3] = Cf[8] + C[0,4] = Cf[9] + C[0,5] = Cf[10] + C[1,2] = Cf[11] + C[1,3] = Cf[12] + C[1,4] = Cf[13] + C[1,5] = Cf[14] + C[2,3] = Cf[15] + C[2,4] = Cf[16] + C[2,5] = Cf[17] + C[3,4] = Cf[18] + C[3,5] = Cf[19] + C[4,5] = Cf[20] + C = np.where(C,C,C.T) + return C + +def compute_delta(): + D = np.zeros((3,3,3,3)) + for a in range(3): + for b in range(3): + for m in range(3): + for n in range(3): + D[a,b,m,n] = (a==m)*(b==n) + (a==n)*(b==m) + d = np.zeros((6,6)) + d[0,0] = D[0,0,0,0] + d[1,1] = D[1,1,1,1] + d[2,2] = D[2,2,2,2] + d[3,3] = D[1,2,1,2] + d[4,4] = D[0,2,0,2] + d[5,5] = D[0,1,0,1] + d[0,1] = D[0,0,1,1] + d[0,2] = D[0,0,2,2] + d[0,3] = D[0,0,1,2] + d[0,4] = D[0,0,0,2] + d[0,5] = D[0,0,0,1] + d[1,2] = D[1,1,2,2] + d[1,3] = D[1,1,1,2] + d[1,4] = D[1,1,0,2] + d[1,5] = D[1,1,0,1] + d[2,3] = D[2,2,1,2] + d[2,4] = D[2,2,0,2] + d[2,5] = D[2,2,0,1] + d[3,4] = D[1,2,0,2] + d[3,5] = D[1,2,0,1] + d[4,5] = D[0,2,0,1] + d = np.where(d,d,d.T) + return d + + +def write_matrix(C, filename): + with open(filename, 'w') as f: + f.write("# Cij Matrix from post process computation\n") + for i in C: + f.write("{:8.5f} {:8.5f} {:8.5f} {:8.5f} {:8.5f} {:8.5f}\n".format( + i[0]*10**-9, i[1]*10**-9, i[2]*10**-9, i[3]*10**-9, i[4]*10**-9, i[5]*10**-9, + ) + ) + return + +def main(): + + N = 500 + vol = 27.047271**3 * 10**-30 # m^3 + T = 60 # K + kb = 1.380649 * 10**-23 # J/K + kbT = T*kb # J + kcalmol2J = 4183.9954/(6.022*10**23) + + born = np.loadtxt('born.out') + stre = np.loadtxt('vir.out') + stre[:, 1:] = -stre[:, 1:]*101325 # -> Pa + try: + mean_born = np.mean(born[:, 1:], axis=0) + except IndexError: + mean_born = born[1:] + + CB = kcalmol2J/vol*reduce_Born(mean_born) # -> J/m^3=Pa + Cs = vol/kbT*np.cov(stre[:,1:].T) + Ct = N*kbT/vol * compute_delta() + C = CB - Cs + Ct + write_matrix(CB, 'born_matrix.out') + write_matrix(Cs, 'stre_matrix.out') + write_matrix(Ct, 'temp_matrix.out') + write_matrix(C, 'full_matrix.out') + C11 = np.mean([C[0,0], C[1,1], C[2,2]]) * 10**-9 + C12 = np.mean([C[0,1], C[0,2], C[1,2]]) * 10**-9 + C44 = np.mean([C[3,3], C[4,4], C[5,5]]) * 10**-9 + eC11 = np.std([C[0,0], C[1,1], C[2,2]]) * 10**-9 + eC12 = np.std([C[0,1], C[0,2], C[1,2]]) * 10**-9 + eC44 = np.std([C[3,3], C[4,4], C[5,5]]) * 10**-9 + print(C*10**-9) + print("C11 = {:f} ± {:f}; C12 = {:f} ± {:f}; C44 = {:f} ± {:f}".format(C11, eC11, C12, eC12, C44, eC44)) + + return + + +if __name__ == "__main__": + try: + main() + except KeyboardInterrupt: + raise SystemExit("User interruption.") + diff --git a/examples/ELASTIC_T/BORN_MATRIX/Argon/Analytical/full_matrix.out b/examples/ELASTIC_T/BORN_MATRIX/Argon/Analytical/full_matrix.out new file mode 100644 index 0000000000..7eb4f6f540 --- /dev/null +++ b/examples/ELASTIC_T/BORN_MATRIX/Argon/Analytical/full_matrix.out @@ -0,0 +1,7 @@ +# Cij Matrix from post process computation + 2.18161 1.13726 1.16596 -0.01607 -0.02637 0.00291 + 1.13726 2.20242 1.16714 0.00386 -0.05820 0.02644 + 1.16596 1.16714 2.24704 -0.00354 -0.00368 0.02714 +-0.01607 0.00386 -0.00354 1.43706 0.00210 0.01003 +-0.02637 -0.05820 -0.00368 0.00210 1.37530 0.01401 + 0.00291 0.02644 0.02714 0.01003 0.01401 1.42403 diff --git a/examples/ELASTIC_T/BORN_MATRIX/Argon/Analytical/in.ljcov b/examples/ELASTIC_T/BORN_MATRIX/Argon/Analytical/in.ljcov new file mode 100644 index 0000000000..255dd65a14 --- /dev/null +++ b/examples/ELASTIC_T/BORN_MATRIX/Argon/Analytical/in.ljcov @@ -0,0 +1,154 @@ +# Analytical calculation +# of Born matrix + +# Note that because of cubic symmetry and central forces, we have: +# C11, pure axial == positive mean value: 1,2,3 +# C44==C23, pure shear == positive mean value, exactly match in pairs: (4,12),(5,8),(6,7) +# C14==C56, shear/axial(normal) == zero mean, exactly match in pairs: (9,21),(14,20),(18,19) +# C15, shear/axial(in-plane) == zero mean: 10,11,13,15,16,17 + +# adjustable parameters + +units real +variable nsteps index 100000 # length of run +variable nthermo index 1000 # thermo output interval +variable nlat equal 5 # size of box +variable T equal 60. # Temperature in K +variable rho equal 5.405 # Lattice spacing in A + +atom_style atomic + +lattice fcc ${rho} +region box block 0 ${nlat} 0 ${nlat} 0 ${nlat} +create_box 1 box +create_atoms 1 box + +mass * 39.948 + +velocity all create ${T} 87287 loop geom +velocity all zero linear + +pair_style lj/cut 12.0 +pair_coeff 1 1 0.238067 3.405 + +neighbor 0.0 bin +neigh_modify every 1 delay 0 check no + +variable vol equal vol +thermo 100 +fix aL all ave/time 1 1 1 v_vol ave running +fix NPT all npt temp $T $T 100 aniso 1. 1. 1000 fixedpoint 0. 0. 0. + +run 20000 + +unfix NPT + +variable newL equal "f_aL^(1./3.)" +change_box all x final 0 ${newL} y final 0. ${newL} z final 0. ${newL} remap units box + +unfix aL + +reset_timestep 0 + +# Conversion variables +variable kb equal 1.38065e-23 # J/K +variable Myvol equal "vol*10^-30" # Volume in m^3 +variable kbt equal "v_kb*v_T" +variable Nat equal atoms +variable Rhokbt equal "v_kbt*v_Nat/v_Myvol" +variable at2Pa equal 101325 +variable kcalmol2J equal "4183.9954/(6.022e23)" +variable C1 equal "v_kcalmol2J/v_Myvol" # Convert Cb from energy to pressure units +variable C2 equal "v_Myvol/v_kbt" # Factor for Cfl terms +variable Pa2GPa equal 1e-9 + +# Born compute giving terms +compute born all born/matrix +# The six virial stress component to compute +compute VIR all pressure NULL virial +variable s1 equal "-c_VIR[1]*v_at2Pa" +variable s2 equal "-c_VIR[2]*v_at2Pa" +variable s3 equal "-c_VIR[3]*v_at2Pa" +variable s6 equal "-c_VIR[4]*v_at2Pa" +variable s5 equal "-c_VIR[5]*v_at2Pa" +variable s4 equal "-c_VIR[6]*v_at2Pa" +variable press equal press + + +# Average of Born term and vector to store stress +# for post processing +fix CB all ave/time 1 ${nthermo} ${nthermo} c_born[*] ave running file born.out overwrite +fix CPR all ave/time 1 1 1 c_VIR[*] file vir.out +fix APR all ave/time 1 1 1 v_press ave running +fix VEC all vector 1 v_s1 v_s2 v_s3 v_s4 v_s5 v_s6 + +thermo ${nthermo} +thermo_style custom step temp press f_APR c_born[1] f_CB[1] c_born[12] f_CB[12] c_born[4] f_CB[4] +thermo_modify line multi + +fix 1 all nvt temp $T $T 100 + +run ${nsteps} + +# Compute vector averages +# Note the indice switch. +# LAMMPS convention is NOT the Voigt notation. +variable aves1 equal "ave(f_VEC[1])" +variable aves2 equal "ave(f_VEC[2])" +variable aves3 equal "ave(f_VEC[3])" +variable aves4 equal "ave(f_VEC[6])" +variable aves5 equal "ave(f_VEC[5])" +variable aves6 equal "ave(f_VEC[4])" + +# Computing the covariance through the - +# is numerically instable. Here we go through the <(s-)^2> +# definition. + +# Computing difference relative to average values +variable ds1 vector "f_VEC[1]-v_aves1" +variable ds2 vector "f_VEC[2]-v_aves2" +variable ds3 vector "f_VEC[3]-v_aves3" +variable ds4 vector "f_VEC[4]-v_aves4" +variable ds5 vector "f_VEC[5]-v_aves5" +variable ds6 vector "f_VEC[6]-v_aves6" + +# Squaring and averaging +variable dds1 vector "v_ds1*v_ds1" +variable dds2 vector "v_ds2*v_ds2" +variable dds3 vector "v_ds3*v_ds3" +variable vars1 equal "ave(v_dds1)" +variable vars2 equal "ave(v_dds2)" +variable vars3 equal "ave(v_dds3)" +variable C11 equal "v_Pa2GPa*(v_C1*f_CB[1] - v_C2*v_vars1 + 2*v_Rhokbt)" +variable C22 equal "v_Pa2GPa*(v_C1*f_CB[2] - v_C2*v_vars2 + 2*v_Rhokbt)" +variable C33 equal "v_Pa2GPa*(v_C1*f_CB[3] - v_C2*v_vars3 + 2*v_Rhokbt)" + +variable dds12 vector "v_ds1*v_ds2" +variable dds13 vector "v_ds1*v_ds3" +variable dds23 vector "v_ds2*v_ds3" +variable vars12 equal "ave(v_dds12)" +variable vars13 equal "ave(v_dds13)" +variable vars23 equal "ave(v_dds23)" +variable C12 equal "v_Pa2GPa*(v_C1*f_CB[7] - v_C2*v_vars12)" +variable C13 equal "v_Pa2GPa*(v_C1*f_CB[8] - v_C2*v_vars13)" +variable C23 equal "v_Pa2GPa*(v_C1*f_CB[12] - v_C2*v_vars23)" + +variable dds4 vector "v_ds4*v_ds4" +variable dds5 vector "v_ds5*v_ds5" +variable dds6 vector "v_ds6*v_ds6" +variable vars4 equal "ave(v_dds4)" +variable vars5 equal "ave(v_dds5)" +variable vars6 equal "ave(v_dds6)" +variable C44 equal "v_Pa2GPa*(v_C1*f_CB[4] - v_C2*v_vars4 + v_Rhokbt)" +variable C55 equal "v_Pa2GPa*(v_C1*f_CB[5] - v_C2*v_vars5 + v_Rhokbt)" +variable C66 equal "v_Pa2GPa*(v_C1*f_CB[6] - v_C2*v_vars6 + v_Rhokbt)" + +variable aC11 equal "(v_C11 + v_C22 + v_C33)/3." +variable aC12 equal "(v_C12 + v_C13 + v_C23)/3." +variable aC44 equal "(v_C44 + v_C55 + v_C66)/3." + +print """ +C11 = ${aC11} +C12 = ${aC12} +C44 = ${aC44} +""" diff --git a/examples/ELASTIC_T/BORN_MATRIX/Argon/Analytical/stre_matrix.out b/examples/ELASTIC_T/BORN_MATRIX/Argon/Analytical/stre_matrix.out new file mode 100644 index 0000000000..94c620b646 --- /dev/null +++ b/examples/ELASTIC_T/BORN_MATRIX/Argon/Analytical/stre_matrix.out @@ -0,0 +1,7 @@ +# Cij Matrix from post process computation + 1.22342 0.73647 0.71011 0.01261 0.02465 -0.00395 + 0.73647 1.20115 0.70711 0.00057 0.05854 -0.02630 + 0.71011 0.70711 1.16055 0.00497 0.00524 -0.02587 + 0.01261 0.00057 0.00497 0.45813 -0.00083 -0.00969 + 0.02465 0.05854 0.00524 -0.00083 0.52170 -0.01747 +-0.00395 -0.02630 -0.02587 -0.00969 -0.01747 0.47064 diff --git a/examples/ELASTIC_T/BORN_MATRIX/Argon/Analytical/temp_matrix.out b/examples/ELASTIC_T/BORN_MATRIX/Argon/Analytical/temp_matrix.out new file mode 100644 index 0000000000..62a734b707 --- /dev/null +++ b/examples/ELASTIC_T/BORN_MATRIX/Argon/Analytical/temp_matrix.out @@ -0,0 +1,7 @@ +# Cij Matrix from post process computation + 0.04187 0.00000 0.00000 0.00000 0.00000 0.00000 + 0.00000 0.04187 0.00000 0.00000 0.00000 0.00000 + 0.00000 0.00000 0.04187 0.00000 0.00000 0.00000 + 0.00000 0.00000 0.00000 0.02093 0.00000 0.00000 + 0.00000 0.00000 0.00000 0.00000 0.02093 0.00000 + 0.00000 0.00000 0.00000 0.00000 0.00000 0.02093 diff --git a/examples/ELASTIC_T/BORN_MATRIX/Argon/Numdiff/born_matrix.out b/examples/ELASTIC_T/BORN_MATRIX/Argon/Numdiff/born_matrix.out new file mode 100644 index 0000000000..01c74114f8 --- /dev/null +++ b/examples/ELASTIC_T/BORN_MATRIX/Argon/Numdiff/born_matrix.out @@ -0,0 +1,7 @@ +# Cij Matrix from post process computation + 3.35855 1.86892 1.87139 0.00233 0.00218 -0.00179 + 1.86892 3.37104 1.87285 0.00112 0.00085 -0.00007 + 1.87139 1.87285 3.37707 -0.00058 0.00038 -0.00057 + 0.00233 0.00112 -0.00058 1.88326 -0.00039 0.00065 + 0.00218 0.00085 0.00038 -0.00039 1.88229 0.00242 +-0.00179 -0.00007 -0.00057 0.00065 0.00242 1.87968 diff --git a/examples/ELASTIC_T/BORN_MATRIX/Argon/Numdiff/compute_born.py b/examples/ELASTIC_T/BORN_MATRIX/Argon/Numdiff/compute_born.py new file mode 100644 index 0000000000..0f3265faeb --- /dev/null +++ b/examples/ELASTIC_T/BORN_MATRIX/Argon/Numdiff/compute_born.py @@ -0,0 +1,118 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +import sys +import numpy as np + +def reduce_Born(Cf): + C = np.zeros((6,6), dtype=np.float64) + C[0,0] = Cf[0] + C[1,1] = Cf[1] + C[2,2] = Cf[2] + C[3,3] = Cf[3] + C[4,4] = Cf[4] + C[5,5] = Cf[5] + C[0,1] = Cf[6] + C[0,2] = Cf[7] + C[0,3] = Cf[8] + C[0,4] = Cf[9] + C[0,5] = Cf[10] + C[1,2] = Cf[11] + C[1,3] = Cf[12] + C[1,4] = Cf[13] + C[1,5] = Cf[14] + C[2,3] = Cf[15] + C[2,4] = Cf[16] + C[2,5] = Cf[17] + C[3,4] = Cf[18] + C[3,5] = Cf[19] + C[4,5] = Cf[20] + C = np.where(C,C,C.T) + return C + +def compute_delta(): + D = np.zeros((3,3,3,3)) + for a in range(3): + for b in range(3): + for m in range(3): + for n in range(3): + D[a,b,m,n] = (a==m)*(b==n) + (a==n)*(b==m) + d = np.zeros((6,6)) + d[0,0] = D[0,0,0,0] + d[1,1] = D[1,1,1,1] + d[2,2] = D[2,2,2,2] + d[3,3] = D[1,2,1,2] + d[4,4] = D[0,2,0,2] + d[5,5] = D[0,1,0,1] + d[0,1] = D[0,0,1,1] + d[0,2] = D[0,0,2,2] + d[0,3] = D[0,0,1,2] + d[0,4] = D[0,0,0,2] + d[0,5] = D[0,0,0,1] + d[1,2] = D[1,1,2,2] + d[1,3] = D[1,1,1,2] + d[1,4] = D[1,1,0,2] + d[1,5] = D[1,1,0,1] + d[2,3] = D[2,2,1,2] + d[2,4] = D[2,2,0,2] + d[2,5] = D[2,2,0,1] + d[3,4] = D[1,2,0,2] + d[3,5] = D[1,2,0,1] + d[4,5] = D[0,2,0,1] + d = np.where(d,d,d.T) + return d + + +def write_matrix(C, filename): + with open(filename, 'w') as f: + f.write("# Cij Matrix from post process computation\n") + for i in C: + f.write("{:8.5f} {:8.5f} {:8.5f} {:8.5f} {:8.5f} {:8.5f}\n".format( + i[0]*10**-9, i[1]*10**-9, i[2]*10**-9, i[3]*10**-9, i[4]*10**-9, i[5]*10**-9, + ) + ) + return + +def main(): + + N = 500 + vol = 27.047271**3 * 10**-30 # m^3 + T = 60 # K + kb = 1.380649 * 10**-23 # J/K + kbT = T*kb # J + kcalmol2J = 4183.9954/(6.022*10**23) + + born = np.loadtxt('born.out') + stre = np.loadtxt('vir.out') + stre[:, 1:] = -stre[:, 1:]*101325 # -> Pa + try: + mean_born = np.mean(born[:, 1:], axis=0) + except IndexError: + mean_born = born[1:] + + CB = kcalmol2J/vol*reduce_Born(mean_born) # -> J/m^3=Pa + Cs = vol/kbT*np.cov(stre[:,1:].T) + Ct = N*kbT/vol * compute_delta() + C = CB - Cs + Ct + write_matrix(CB, 'born_matrix.out') + write_matrix(Cs, 'stre_matrix.out') + write_matrix(Ct, 'temp_matrix.out') + write_matrix(C, 'full_matrix.out') + C11 = np.mean([C[0,0], C[1,1], C[2,2]]) * 10**-9 + C12 = np.mean([C[0,1], C[0,2], C[1,2]]) * 10**-9 + C44 = np.mean([C[3,3], C[4,4], C[5,5]]) * 10**-9 + eC11 = np.std([C[0,0], C[1,1], C[2,2]]) * 10**-9 + eC12 = np.std([C[0,1], C[0,2], C[1,2]]) * 10**-9 + eC44 = np.std([C[3,3], C[4,4], C[5,5]]) * 10**-9 + print(C*10**-9) + print("C11 = {:f} ± {:f}; C12 = {:f} ± {:f}; C44 = {:f} ± {:f}".format(C11, eC11, C12, eC12, C44, eC44)) + + return + + +if __name__ == "__main__": + try: + main() + except KeyboardInterrupt: + raise SystemExit("User interruption.") + diff --git a/examples/ELASTIC_T/BORN_MATRIX/Argon/Numdiff/full_matrix.out b/examples/ELASTIC_T/BORN_MATRIX/Argon/Numdiff/full_matrix.out new file mode 100644 index 0000000000..93a2c45e42 --- /dev/null +++ b/examples/ELASTIC_T/BORN_MATRIX/Argon/Numdiff/full_matrix.out @@ -0,0 +1,7 @@ +# Cij Matrix from post process computation + 2.29922 1.17439 1.20854 -0.01653 -0.01684 0.01188 + 1.17439 2.20673 1.21718 -0.00781 -0.00753 0.00867 + 1.20854 1.21718 2.30804 0.01535 -0.01596 0.00426 +-0.01653 -0.00781 0.01535 1.47647 -0.01355 -0.01601 +-0.01684 -0.00753 -0.01596 -0.01355 1.37905 0.01975 + 0.01188 0.00867 0.00426 -0.01601 0.01975 1.40170 diff --git a/examples/ELASTIC_T/BORN_MATRIX/Argon/Numdiff/in.ljcov b/examples/ELASTIC_T/BORN_MATRIX/Argon/Numdiff/in.ljcov new file mode 100644 index 0000000000..6db9cfaa66 --- /dev/null +++ b/examples/ELASTIC_T/BORN_MATRIX/Argon/Numdiff/in.ljcov @@ -0,0 +1,154 @@ +# Numerical difference calculation +# of Born matrix + +# Note that because of cubic symmetry and central forces, we have: +# C11, pure axial == positive mean value: 1,2,3 +# C44==C23, pure shear == positive mean value, exactly match in pairs: (4,12),(5,8),(6,7) +# C14==C56, shear/axial(normal) == zero mean, exactly match in pairs: (9,21),(14,20),(18,19) +# C15, shear/axial(in-plane) == zero mean: 10,11,13,15,16,17 + +# adjustable parameters + +units real +variable nsteps index 100000 # length of run +variable nthermo index 1000 # thermo output interval +variable nlat equal 5 # size of box +variable T equal 60. # Temperature in K +variable rho equal 5.405 # Lattice spacing in A + +atom_style atomic + +lattice fcc ${rho} +region box block 0 ${nlat} 0 ${nlat} 0 ${nlat} +create_box 1 box +create_atoms 1 box + +mass * 39.948 + +velocity all create ${T} 87287 loop geom +velocity all zero linear + +pair_style lj/cut 12.0 +pair_coeff 1 1 0.238067 3.405 + +neighbor 0.0 bin +neigh_modify every 1 delay 0 check no + +variable vol equal vol +thermo 100 +fix aL all ave/time 1 1 1 v_vol ave running +fix NPT all npt temp $T $T 100 aniso 1. 1. 1000 fixedpoint 0. 0. 0. + +run 20000 + +unfix NPT + +variable newL equal "f_aL^(1./3.)" +change_box all x final 0 ${newL} y final 0. ${newL} z final 0. ${newL} remap units box + +unfix aL + +reset_timestep 0 + +# Conversion variables +variable kb equal 1.38065e-23 # J/K +variable Myvol equal "vol*10^-30" # Volume in m^3 +variable kbt equal "v_kb*v_T" +variable Nat equal atoms +variable Rhokbt equal "v_kbt*v_Nat/v_Myvol" +variable at2Pa equal 101325 +variable kcalmol2J equal "4183.9954/(6.022e23)" +variable C1 equal "v_kcalmol2J/v_Myvol" # Convert Cb from energy to pressure units +variable C2 equal "v_Myvol/v_kbt" # Factor for Cfl terms +variable Pa2GPa equal 1e-9 + +# Born compute giving terms +# The six virial stress component to compute +compute VIR all pressure NULL virial +compute born all born/matrix numdiff 1e-6 VIR +variable s1 equal "-c_VIR[1]*v_at2Pa" +variable s2 equal "-c_VIR[2]*v_at2Pa" +variable s3 equal "-c_VIR[3]*v_at2Pa" +variable s6 equal "-c_VIR[4]*v_at2Pa" +variable s5 equal "-c_VIR[5]*v_at2Pa" +variable s4 equal "-c_VIR[6]*v_at2Pa" +variable press equal press + + +# Average of Born term and vector to store stress +# for post processing +fix CB all ave/time 1 ${nthermo} ${nthermo} c_born[*] ave running file born.out overwrite +fix CPR all ave/time 1 1 1 c_VIR[*] file vir.out +fix APR all ave/time 1 1 1 v_press ave running +fix VEC all vector 1 v_s1 v_s2 v_s3 v_s4 v_s5 v_s6 + +thermo ${nthermo} +thermo_style custom step temp press f_APR c_born[1] f_CB[1] c_born[12] f_CB[12] c_born[4] f_CB[4] +thermo_modify line multi + +fix 1 all nvt temp $T $T 100 + +run ${nsteps} + +# Compute vector averages +# Note the indice switch. +# LAMMPS convention is NOT the Voigt notation. +variable aves1 equal "ave(f_VEC[1])" +variable aves2 equal "ave(f_VEC[2])" +variable aves3 equal "ave(f_VEC[3])" +variable aves4 equal "ave(f_VEC[6])" +variable aves5 equal "ave(f_VEC[5])" +variable aves6 equal "ave(f_VEC[4])" + +# Computing the covariance through the - +# is numerically instable. Here we go through the <(s-)^2> +# definition. + +# Computing difference relative to average values +variable ds1 vector "f_VEC[1]-v_aves1" +variable ds2 vector "f_VEC[2]-v_aves2" +variable ds3 vector "f_VEC[3]-v_aves3" +variable ds4 vector "f_VEC[4]-v_aves4" +variable ds5 vector "f_VEC[5]-v_aves5" +variable ds6 vector "f_VEC[6]-v_aves6" + +# Squaring and averaging +variable dds1 vector "v_ds1*v_ds1" +variable dds2 vector "v_ds2*v_ds2" +variable dds3 vector "v_ds3*v_ds3" +variable vars1 equal "ave(v_dds1)" +variable vars2 equal "ave(v_dds2)" +variable vars3 equal "ave(v_dds3)" +variable C11 equal "v_Pa2GPa*(v_C1*f_CB[1] - v_C2*v_vars1 + 2*v_Rhokbt)" +variable C22 equal "v_Pa2GPa*(v_C1*f_CB[2] - v_C2*v_vars2 + 2*v_Rhokbt)" +variable C33 equal "v_Pa2GPa*(v_C1*f_CB[3] - v_C2*v_vars3 + 2*v_Rhokbt)" + +variable dds12 vector "v_ds1*v_ds2" +variable dds13 vector "v_ds1*v_ds3" +variable dds23 vector "v_ds2*v_ds3" +variable vars12 equal "ave(v_dds12)" +variable vars13 equal "ave(v_dds13)" +variable vars23 equal "ave(v_dds23)" +variable C12 equal "v_Pa2GPa*(v_C1*f_CB[7] - v_C2*v_vars12)" +variable C13 equal "v_Pa2GPa*(v_C1*f_CB[8] - v_C2*v_vars13)" +variable C23 equal "v_Pa2GPa*(v_C1*f_CB[12] - v_C2*v_vars23)" + +variable dds4 vector "v_ds4*v_ds4" +variable dds5 vector "v_ds5*v_ds5" +variable dds6 vector "v_ds6*v_ds6" +variable vars4 equal "ave(v_dds4)" +variable vars5 equal "ave(v_dds5)" +variable vars6 equal "ave(v_dds6)" +variable C44 equal "v_Pa2GPa*(v_C1*f_CB[4] - v_C2*v_vars4 + v_Rhokbt)" +variable C55 equal "v_Pa2GPa*(v_C1*f_CB[5] - v_C2*v_vars5 + v_Rhokbt)" +variable C66 equal "v_Pa2GPa*(v_C1*f_CB[6] - v_C2*v_vars6 + v_Rhokbt)" + +variable aC11 equal "(v_C11 + v_C22 + v_C33)/3." +variable aC12 equal "(v_C12 + v_C13 + v_C23)/3." +variable aC44 equal "(v_C44 + v_C55 + v_C66)/3." + +print """ +C11 = ${aC11} +C12 = ${aC12} +C44 = ${aC44} +""" diff --git a/examples/ELASTIC_T/BORN_MATRIX/Argon/Numdiff/stre_matrix.out b/examples/ELASTIC_T/BORN_MATRIX/Argon/Numdiff/stre_matrix.out new file mode 100644 index 0000000000..c6d64581a2 --- /dev/null +++ b/examples/ELASTIC_T/BORN_MATRIX/Argon/Numdiff/stre_matrix.out @@ -0,0 +1,7 @@ +# Cij Matrix from post process computation + 1.10120 0.69454 0.66285 0.01885 0.01902 -0.01367 + 0.69454 1.20617 0.65567 0.00893 0.00839 -0.00873 + 0.66285 0.65567 1.11090 -0.01593 0.01634 -0.00483 + 0.01885 0.00893 -0.01593 0.42772 0.01316 0.01666 + 0.01902 0.00839 0.01634 0.01316 0.52416 -0.01733 +-0.01367 -0.00873 -0.00483 0.01666 -0.01733 0.49891 diff --git a/examples/ELASTIC_T/BORN_MATRIX/Argon/Numdiff/temp_matrix.out b/examples/ELASTIC_T/BORN_MATRIX/Argon/Numdiff/temp_matrix.out new file mode 100644 index 0000000000..62a734b707 --- /dev/null +++ b/examples/ELASTIC_T/BORN_MATRIX/Argon/Numdiff/temp_matrix.out @@ -0,0 +1,7 @@ +# Cij Matrix from post process computation + 0.04187 0.00000 0.00000 0.00000 0.00000 0.00000 + 0.00000 0.04187 0.00000 0.00000 0.00000 0.00000 + 0.00000 0.00000 0.04187 0.00000 0.00000 0.00000 + 0.00000 0.00000 0.00000 0.02093 0.00000 0.00000 + 0.00000 0.00000 0.00000 0.00000 0.02093 0.00000 + 0.00000 0.00000 0.00000 0.00000 0.00000 0.02093 diff --git a/examples/ELASTIC_T/BORN_MATRIX/Argon/README.md b/examples/ELASTIC_T/BORN_MATRIX/Argon/README.md new file mode 100644 index 0000000000..c0a11af16f --- /dev/null +++ b/examples/ELASTIC_T/BORN_MATRIX/Argon/README.md @@ -0,0 +1,13 @@ +This repository is a test case for the compute born/matrix. It provides short +scripts creating argon fcc crystal and computing the Born term using the two +methods described in the documentation. + +In the __Analytical__ directory the terms are computed using the analytical +derivation of the Born term for the lj/cut pair style. + +In the __Numdiff__ directory, the Born term is evaluated through small +numerical differences of the stress tensor. This method can be used with any +interaction potential. + +Both script show examples on how to compute the full Cij elastic stiffness +tensor in LAMMPS. diff --git a/examples/ELASTIC_T/BORN_MATRIX/Silicon/Si.sw b/examples/ELASTIC_T/BORN_MATRIX/Silicon/Si.sw new file mode 120000 index 0000000000..e575921334 --- /dev/null +++ b/examples/ELASTIC_T/BORN_MATRIX/Silicon/Si.sw @@ -0,0 +1 @@ +../../../../potentials/Si.sw \ No newline at end of file diff --git a/examples/ELASTIC_T/BORN_MATRIX/Silicon/final_output.in b/examples/ELASTIC_T/BORN_MATRIX/Silicon/final_output.in new file mode 100644 index 0000000000..e0d7777b4b --- /dev/null +++ b/examples/ELASTIC_T/BORN_MATRIX/Silicon/final_output.in @@ -0,0 +1,68 @@ +# Average moduli for cubic crystals + +variable C11cubic equal (${C11}+${C22}+${C33})/3.0 +variable C12cubic equal (${C12}+${C13}+${C23})/3.0 +variable C44cubic equal (${C44}+${C55}+${C66})/3.0 + +variable bulkmodulus equal (${C11cubic}+2*${C12cubic})/3.0 +variable shearmodulus1 equal ${C44cubic} +variable shearmodulus2 equal (${C11cubic}-${C12cubic})/2.0 +variable poissonratio equal 1.0/(1.0+${C11cubic}/${C12cubic}) + +# For Stillinger-Weber silicon, the analytical results +# are known to be (E. R. Cowley, 1988): +# C11 = 151.4 GPa +# C12 = 76.4 GPa +# C44 = 56.4 GPa + +#print "=========================================" +#print "Components of the Elastic Constant Tensor" +#print "=========================================" + +print "Elastic Constant C11 = ${C11} ${cunits}" +print "Elastic Constant C22 = ${C22} ${cunits}" +print "Elastic Constant C33 = ${C33} ${cunits}" + +print "Elastic Constant C12 = ${C12} ${cunits}" +print "Elastic Constant C13 = ${C13} ${cunits}" +print "Elastic Constant C23 = ${C23} ${cunits}" + +print "Elastic Constant C44 = ${C44} ${cunits}" +print "Elastic Constant C55 = ${C55} ${cunits}" +print "Elastic Constant C66 = ${C66} ${cunits}" + +print "Elastic Constant C14 = ${C14} ${cunits}" +print "Elastic Constant C15 = ${C15} ${cunits}" +print "Elastic Constant C16 = ${C16} ${cunits}" + +print "Elastic Constant C24 = ${C24} ${cunits}" +print "Elastic Constant C25 = ${C25} ${cunits}" +print "Elastic Constant C26 = ${C26} ${cunits}" + +print "Elastic Constant C34 = ${C34} ${cunits}" +print "Elastic Constant C35 = ${C35} ${cunits}" +print "Elastic Constant C36 = ${C36} ${cunits}" + +print "Elastic Constant C45 = ${C45} ${cunits}" +print "Elastic Constant C46 = ${C46} ${cunits}" +print "Elastic Constant C56 = ${C56} ${cunits}" + +print "=========================================" +print "Average properties for a cubic crystal" +print "=========================================" + +print "Bulk Modulus = ${bulkmodulus} ${cunits}" +print "Shear Modulus 1 = ${shearmodulus1} ${cunits}" +print "Shear Modulus 2 = ${shearmodulus2} ${cunits}" +print "Poisson Ratio = ${poissonratio}" + +# summarize sampling protocol + +variable tmp equal atoms +print "Number of atoms = ${tmp}" +print "Stress sampling interval = ${nevery}" +variable tmp equal ${nrun}/${nevery} +print "Stress sample count = ${tmp}" +print "Born sampling interval = ${neveryborn}" +variable tmp equal ${nrun}/${neveryborn} +print "Born sample count = ${tmp}" diff --git a/examples/ELASTIC_T/BORN_MATRIX/Silicon/in.elastic b/examples/ELASTIC_T/BORN_MATRIX/Silicon/in.elastic new file mode 100644 index 0000000000..cbd0aa208c --- /dev/null +++ b/examples/ELASTIC_T/BORN_MATRIX/Silicon/in.elastic @@ -0,0 +1,25 @@ +# Compute elastic constant tensor for a crystal at finite temperature +# These settings replicate the 1477~K benchmark of +# Kluge, Ray, and Rahman (1986) that is Ref.[15] in: +# Y. Zhen, C. Chu, Computer Physics Communications 183(2012) 261-265 + +# here: Y. Zhen, C. Chu, Computer Physics Communications 183(2012) 261-265 + +include init.in + +# Compute initial state + +include potential.in +thermo_style custom step temp pe press density +run ${nequil} + +# Run dynamics + +include potential.in +include output.in + +run ${nrun} + +# Output final values + +include final_output.in diff --git a/examples/ELASTIC_T/BORN_MATRIX/Silicon/init.in b/examples/ELASTIC_T/BORN_MATRIX/Silicon/init.in new file mode 100644 index 0000000000..09aa85a489 --- /dev/null +++ b/examples/ELASTIC_T/BORN_MATRIX/Silicon/init.in @@ -0,0 +1,59 @@ +# NOTE: This script can be modified for different atomic structures, +# units, etc. See in.elastic for more info. +# + +# Define MD parameters +# These can be modified by the user +# These settings replicate the 1477~K benchmark of +# Kluge, Ray, and Rahman (1986) that is Ref.[15] in: +# Y. Zhen, C. Chu, Computer Physics Communications 183(2012) 261-265 + +# select temperature and pressure (lattice constant) + +variable temp index 1477.0 # temperature of initial sample +variable a index 5.457 # lattice constant + +# select sampling parameters, important for speed/convergence + +variable nthermo index 1500 # interval for thermo output +variable nevery index 10 # stress sampling interval +variable neveryborn index 100 # Born sampling interval +variable timestep index 0.000766 # timestep +variable nlat index 3 # number of lattice unit cells + +# other settings + +variable mass1 index 28.06 # mass +variable tdamp index 0.01 # time constant for thermostat +variable seed index 123457 # seed for thermostat +variable thermostat index 1 # 0 if NVE, 1 if NVT +variable delta index 1.0e-6 # Born numdiff strain magnitude + +# hard-coded rules-of-thumb for run length, etc. + +variable nfreq equal ${nthermo} # interval for averaging output +variable nrepeat equal floor(${nfreq}/${nevery}) # number of samples +variable nrepeatborn equal floor(${nfreq}/${neveryborn}) # number of samples +variable nequil equal 10*${nthermo} # length of equilibration run +variable nrun equal 100*${nthermo} # length of equilibrated run + +# generate the box and atom positions using a diamond lattice + +units metal + +boundary p p p + +# this generates a standard 8-atom cubic cell + +lattice diamond $a +region box prism 0 1 0 1 0 1 0 0 0 + +# this generates a 2-atom triclinic cell +#include tri.in + +create_box 1 box +create_atoms 1 box +mass 1 ${mass1} +replicate ${nlat} ${nlat} ${nlat} +velocity all create ${temp} 87287 + diff --git a/examples/ELASTIC_T/BORN_MATRIX/Silicon/output.in b/examples/ELASTIC_T/BORN_MATRIX/Silicon/output.in new file mode 100644 index 0000000000..d6a89cff6a --- /dev/null +++ b/examples/ELASTIC_T/BORN_MATRIX/Silicon/output.in @@ -0,0 +1,140 @@ +# Setup output + +# Stress fluctuation term F + +compute stress all pressure thermo_temp +variable s1 equal c_stress[1] +variable s2 equal c_stress[2] +variable s3 equal c_stress[3] +variable s4 equal c_stress[6] +variable s5 equal c_stress[5] +variable s6 equal c_stress[4] + +variable s11 equal v_s1*v_s1 +variable s22 equal v_s2*v_s2 +variable s33 equal v_s3*v_s3 +variable s44 equal v_s4*v_s4 +variable s55 equal v_s5*v_s5 +variable s66 equal v_s6*v_s6 +variable s33 equal v_s3*v_s3 +variable s12 equal v_s1*v_s2 +variable s13 equal v_s1*v_s3 +variable s14 equal v_s1*v_s4 +variable s15 equal v_s1*v_s5 +variable s16 equal v_s1*v_s6 +variable s23 equal v_s2*v_s3 +variable s24 equal v_s2*v_s4 +variable s25 equal v_s2*v_s5 +variable s26 equal v_s2*v_s6 +variable s34 equal v_s3*v_s4 +variable s35 equal v_s3*v_s5 +variable s36 equal v_s3*v_s6 +variable s45 equal v_s4*v_s5 +variable s46 equal v_s4*v_s6 +variable s56 equal v_s5*v_s6 + +variable mytemp equal temp +variable mypress equal press +variable mype equal pe/atoms +fix avt all ave/time ${nevery} ${nrepeat} ${nfreq} v_mytemp ave running +fix avp all ave/time ${nevery} ${nrepeat} ${nfreq} v_mypress ave running +fix avpe all ave/time ${nevery} ${nrepeat} ${nfreq} v_mype ave running +fix avs all ave/time ${nevery} ${nrepeat} ${nfreq} v_s1 v_s2 v_s3 v_s4 v_s5 v_s6 ave running +fix avssq all ave/time ${nevery} ${nrepeat} ${nfreq} & +v_s11 v_s22 v_s33 v_s44 v_s55 v_s66 & +v_s12 v_s13 v_s14 v_s15 v_s16 & +v_s23 v_s24 v_s25 v_s26 & +v_s34 v_s35 v_s36 & +v_s45 v_s46 & +v_s56 & +ave running + +# bar to GPa +variable pconv equal 1.0e5/1.0e9 +variable cunits index GPa +# metal unit constants from LAMMPS +# force->nktv2p = 1.6021765e6; +# force->boltz = 8.617343e-5; +variable boltz equal 8.617343e-5 +variable nktv2p equal 1.6021765e6 +variable vkt equal vol/(${boltz}*${temp})/${nktv2p} +variable ffac equal ${pconv}*${vkt} + +variable F11 equal -(f_avssq[1]-f_avs[1]*f_avs[1])*${ffac} +variable F22 equal -(f_avssq[2]-f_avs[2]*f_avs[2])*${ffac} +variable F33 equal -(f_avssq[3]-f_avs[3]*f_avs[3])*${ffac} +variable F44 equal -(f_avssq[4]-f_avs[4]*f_avs[4])*${ffac} +variable F55 equal -(f_avssq[5]-f_avs[5]*f_avs[5])*${ffac} +variable F66 equal -(f_avssq[6]-f_avs[6]*f_avs[6])*${ffac} + +variable F12 equal -(f_avssq[7]-f_avs[1]*f_avs[2])*${ffac} +variable F13 equal -(f_avssq[8]-f_avs[1]*f_avs[3])*${ffac} +variable F14 equal -(f_avssq[9]-f_avs[1]*f_avs[4])*${ffac} +variable F15 equal -(f_avssq[10]-f_avs[1]*f_avs[5])*${ffac} +variable F16 equal -(f_avssq[11]-f_avs[1]*f_avs[6])*${ffac} + +variable F23 equal -(f_avssq[12]-f_avs[2]*f_avs[3])*${ffac} +variable F24 equal -(f_avssq[13]-f_avs[2]*f_avs[4])*${ffac} +variable F25 equal -(f_avssq[14]-f_avs[2]*f_avs[5])*${ffac} +variable F26 equal -(f_avssq[15]-f_avs[2]*f_avs[6])*${ffac} + +variable F34 equal -(f_avssq[16]-f_avs[3]*f_avs[4])*${ffac} +variable F35 equal -(f_avssq[17]-f_avs[3]*f_avs[5])*${ffac} +variable F36 equal -(f_avssq[18]-f_avs[3]*f_avs[6])*${ffac} + +variable F45 equal -(f_avssq[19]-f_avs[4]*f_avs[5])*${ffac} +variable F46 equal -(f_avssq[20]-f_avs[4]*f_avs[6])*${ffac} + +variable F56 equal -(f_avssq[21]-f_avs[5]*f_avs[6])*${ffac} + +# Born term + +compute virial all pressure NULL virial +compute born all born/matrix numdiff ${delta} virial +fix avborn all ave/time ${neveryborn} ${nrepeatborn} ${nfreq} c_born[*] ave running + +variable bfac equal ${pconv}*${nktv2p}/vol +variable B vector f_avborn*${bfac} + +# Kinetic term + +variable kfac equal ${pconv}*${nktv2p}*atoms*${boltz}*${temp}/vol +variable K11 equal 4.0*${kfac} +variable K22 equal 4.0*${kfac} +variable K33 equal 4.0*${kfac} +variable K44 equal 2.0*${kfac} +variable K55 equal 2.0*${kfac} +variable K66 equal 2.0*${kfac} + +# Add F, K, and B together + +variable C11 equal v_F11+v_B[1]+v_K11 +variable C22 equal v_F22+v_B[2]+v_K22 +variable C33 equal v_F33+v_B[3]+v_K33 +variable C44 equal v_F44+v_B[4]+v_K44 +variable C55 equal v_F55+v_B[5]+v_K55 +variable C66 equal v_F66+v_B[6]+v_K66 + +variable C12 equal v_F12+v_B[7] +variable C13 equal v_F13+v_B[8] +variable C14 equal v_F14+v_B[9] +variable C15 equal v_F15+v_B[10] +variable C16 equal v_F16+v_B[11] + +variable C23 equal v_F23+v_B[12] +variable C24 equal v_F24+v_B[13] +variable C25 equal v_F25+v_B[14] +variable C26 equal v_F26+v_B[15] + +variable C34 equal v_F34+v_B[16] +variable C35 equal v_F35+v_B[17] +variable C36 equal v_F36+v_B[18] + +variable C45 equal v_F45+v_B[19] +variable C46 equal v_F46+v_B[20] + +variable C56 equal v_F56+v_B[21] + +thermo ${nthermo} +thermo_style custom step temp pe press density f_avt f_avp f_avpe v_F11 v_F22 v_F33 v_F44 v_F55 v_F66 v_F12 v_F13 v_F23 v_B[1] v_B[2] v_B[3] v_B[4] v_B[5] v_B[6] v_B[7] v_B[8] v_B[12] +thermo_modify norm no diff --git a/examples/ELASTIC_T/BORN_MATRIX/Silicon/potential.in b/examples/ELASTIC_T/BORN_MATRIX/Silicon/potential.in new file mode 100644 index 0000000000..43b0fa61c4 --- /dev/null +++ b/examples/ELASTIC_T/BORN_MATRIX/Silicon/potential.in @@ -0,0 +1,21 @@ +# NOTE: This script can be modified for different pair styles +# See in.elastic for more info. + +reset_timestep 0 + +# Choose potential +pair_style sw +pair_coeff * * Si.sw Si + +# Setup neighbor style +neighbor 1.0 nsq +neigh_modify once no every 1 delay 0 check yes + +# Setup MD + +timestep ${timestep} +fix 4 all nve +if "${thermostat} == 1" then & + "fix 5 all langevin ${temp} ${temp} ${tdamp} ${seed}" + + diff --git a/examples/ELASTIC_T/BORN_MATRIX/Silicon/tri.in b/examples/ELASTIC_T/BORN_MATRIX/Silicon/tri.in new file mode 100644 index 0000000000..20bfca16ec --- /dev/null +++ b/examples/ELASTIC_T/BORN_MATRIX/Silicon/tri.in @@ -0,0 +1,26 @@ +# this generates a 2-atom triclinic cell +# due to rotation on to x-axis, +# elastic constant analysis is not working yet + +# unit lattice vectors are +# a1 = (1 0 0) +# a2 = (1/2 sqrt3/2 0) +# a3 = (1/2 1/(2sqrt3) sqrt2/sqrt3) + +variable a1x equal 1 +variable a2x equal 1/2 +variable a2y equal sqrt(3)/2 +variable a3x equal 1/2 +variable a3y equal 1/(2*sqrt(3)) +variable a3z equal sqrt(2/3) +variable l equal $a/sqrt(2) + +lattice custom ${l} & + a1 ${a1x} 0 0 & + a2 ${a2x} ${a2y} 0.0 & + a3 ${a3x} ${a3y} ${a3z} & + basis 0 0 0 & + basis 0.25 0.25 0.25 & + spacing 1 1 1 + +region box prism 0 ${a1x} 0 ${a2y} 0 ${a3z} ${a2x} ${a3x} ${a3y} diff --git a/examples/ELASTIC_T/DEFORMATION/Silicon/Si.sw b/examples/ELASTIC_T/DEFORMATION/Silicon/Si.sw new file mode 120000 index 0000000000..e575921334 --- /dev/null +++ b/examples/ELASTIC_T/DEFORMATION/Silicon/Si.sw @@ -0,0 +1 @@ +../../../../potentials/Si.sw \ No newline at end of file diff --git a/examples/ELASTIC_T/displace.mod b/examples/ELASTIC_T/DEFORMATION/Silicon/displace.mod similarity index 100% rename from examples/ELASTIC_T/displace.mod rename to examples/ELASTIC_T/DEFORMATION/Silicon/displace.mod diff --git a/examples/ELASTIC_T/in.elastic b/examples/ELASTIC_T/DEFORMATION/Silicon/in.elastic similarity index 100% rename from examples/ELASTIC_T/in.elastic rename to examples/ELASTIC_T/DEFORMATION/Silicon/in.elastic diff --git a/examples/ELASTIC_T/init.mod b/examples/ELASTIC_T/DEFORMATION/Silicon/init.mod similarity index 100% rename from examples/ELASTIC_T/init.mod rename to examples/ELASTIC_T/DEFORMATION/Silicon/init.mod diff --git a/examples/ELASTIC_T/potential.mod b/examples/ELASTIC_T/DEFORMATION/Silicon/potential.mod similarity index 93% rename from examples/ELASTIC_T/potential.mod rename to examples/ELASTIC_T/DEFORMATION/Silicon/potential.mod index d4b7cc7158..7cd996d5d5 100644 --- a/examples/ELASTIC_T/potential.mod +++ b/examples/ELASTIC_T/DEFORMATION/Silicon/potential.mod @@ -2,7 +2,7 @@ # See in.elastic for more info. # we must undefine any fix ave/* fix before using reset_timestep -if "$(is_defined(fix,avp)" then "unfix avp" +if "$(is_defined(fix,avp))" then "unfix avp" reset_timestep 0 # Choose potential diff --git a/examples/ELASTIC_T/Si.sw b/examples/ELASTIC_T/Si.sw deleted file mode 100644 index db4be100ef..0000000000 --- a/examples/ELASTIC_T/Si.sw +++ /dev/null @@ -1,18 +0,0 @@ -# DATE: 2007-06-11 CONTRIBUTOR: Aidan Thompson, athomps@sandia.gov CITATION: Stillinger and Weber, Phys Rev B, 31, 5262, (1985) -# Stillinger-Weber parameters for various elements and mixtures -# multiple entries can be added to this file, LAMMPS reads the ones it needs -# these entries are in LAMMPS "metal" units: -# epsilon = eV; sigma = Angstroms -# other quantities are unitless - -# format of a single entry (one or more lines): -# element 1, element 2, element 3, -# epsilon, sigma, a, lambda, gamma, costheta0, A, B, p, q, tol - -# Here are the original parameters in metal units, for Silicon from: -# -# Stillinger and Weber, Phys. Rev. B, v. 31, p. 5262, (1985) -# - -Si Si Si 2.1683 2.0951 1.80 21.0 1.20 -0.333333333333 - 7.049556277 0.6022245584 4.0 0.0 0.0 diff --git a/examples/README b/examples/README index f82fe3b929..8d0d7cf1c9 100644 --- a/examples/README +++ b/examples/README @@ -94,7 +94,7 @@ msst: MSST shock dynamics nb3b: use of nonbonded 3-body harmonic pair style neb: nudged elastic band (NEB) calculation for barrier finding nemd: non-equilibrium MD of 2d sheared system -numdiff: numerical difference computation of forces and virial +numdiff: numerical difference computation of forces, virial, and Born matrix obstacle: flow around two voids in a 2d channel peptide: dynamics of a small solvated peptide chain (5-mer) peri: Peridynamic model of cylinder impacted by indenter @@ -153,12 +153,19 @@ either by itself or in tandem with another code or library. See the COUPLE/README file to get started. The ELASTIC directory has an example script for computing elastic -constants at zero temperature, using an Si example. See the +stiffness tensor (elastic constants) +at zero temperature, using an Si example. See the ELASTIC/in.elastic file for more info. -The ELASTIC_T directory has an example script for computing elastic -constants at finite temperature, using an Si example. See the -ELASTIC_T/in.elastic file for more info. +The ELASTIC_T directory has example scripts for the computing elastic +stiffness tensor at finite temperature. Two different methods are +demonstrated. DEFORMATION estimates the change in the average +stress tensor between multiple simulations +in which small finite deformations are made to the simulation cell. +BORN_MATRIX runs a single simulation in which the Born matrix and stress +fluctuations are averaged. The second method +is newer in LAMMPS and is generally more efficient and +more reliable. The HEAT directory has example scripts for heat exchange algorithms (e.g. used for establishing a thermal gradient), using two different diff --git a/examples/numdiff/README.md b/examples/numdiff/README.md index d88f66d65c..10727fe0e3 100644 --- a/examples/numdiff/README.md +++ b/examples/numdiff/README.md @@ -1,8 +1,11 @@ -# LAMMPS FIX NUMDIFF EXAMPLE +# LAMMPS NUMDIFF EXAMPLES FOR FORCES, VIRIAL, and BORN MATRIX -## Numerical Difference Fix +## Numerical Difference Fixes and Computes -This directory contains the ingredients to run an NVE simulation using the numerical difference fix and calculate error in forces. +This directory contains the input script for an NVE simulation with +fix numdiff, fix numdiff/virial, and compute born/matrix numdiff. +In each cases, results are compared to exact analytic expressions +and the small relative differences are reported. Example: ``` @@ -10,4 +13,4 @@ NP=4 #number of processors mpirun -np $NP lmp_mpi -in.numdiff ``` -## Required LAMMPS packages: MOLECULE package +## Required LAMMPS packages: EXTRA-FIX, EXTRA-COMPUTE diff --git a/examples/numdiff/in.numdiff b/examples/numdiff/in.numdiff index cd3d44a9e7..363570da87 100644 --- a/examples/numdiff/in.numdiff +++ b/examples/numdiff/in.numdiff @@ -1,5 +1,5 @@ # Numerical difference calculation -# of error in forces and virial stress +# of error in forces, virial stress, and Born matrix # adjustable parameters @@ -8,8 +8,10 @@ variable nthermo index 10 # thermo output interval variable ndump index 500 # dump output interval variable nlat index 3 # size of box variable fdelta index 1.0e-4 # displacement size -variable vdelta index 1.0e-6 # strain size +variable vdelta index 1.0e-6 # strain size for numdiff/virial +variable bdelta index 1.0e-8 # strain size for numdiff Born matrix variable temp index 10.0 # temperature +variable nugget equal 1.0e-6 # regularization for relerr units metal atom_style atomic @@ -23,8 +25,8 @@ mass 1 39.903 velocity all create ${temp} 2357 mom yes dist gaussian -pair_style lj/cubic -pair_coeff * * 0.0102701 3.42 +pair_style lj/cut 5.0 +pair_coeff 1 1 0.0102701 3.42 neighbor 0.0 bin neigh_modify every 1 delay 0 check no @@ -42,7 +44,7 @@ variable ferrsq atom v_ferrx^2+v_ferry^2+v_ferrz^2 compute faverrsq all reduce ave v_ferrsq variable fsq atom fx^2+fy^2+fz^2 compute favsq all reduce ave v_fsq -variable frelerr equal sqrt(c_faverrsq/c_favsq) +variable frelerr equal sqrt(c_faverrsq/(c_favsq+${nugget})) dump errors all custom ${ndump} force_error.dump v_ferrx v_ferry v_ferrz # define numerical virial stress tensor calculation @@ -67,8 +69,59 @@ variable vsq equal "c_myvirial[1]^2 + & c_myvirial[4]^2 + & c_myvirial[5]^2 + & c_myvirial[6]^2" -variable vrelerr equal sqrt(v_verrsq/v_vsq) +variable vrelerr equal sqrt(v_verrsq/(v_vsq+${nugget})) -thermo_style custom step temp pe etotal press v_frelerr v_vrelerr +# define numerical Born matrix calculation + +compute bornnum all born/matrix numdiff ${bdelta} myvirial +compute born all born/matrix +variable berr vector c_bornnum-c_born +variable berrsq equal "v_berr[1]^2 + & + v_berr[2]^2 + & + v_berr[3]^2 + & + v_berr[4]^2 + & + v_berr[5]^2 + & + v_berr[6]^2 + & + v_berr[7]^2 + & + v_berr[8]^2 + & + v_berr[9]^2 + & + v_berr[10]^2 + & + v_berr[11]^2 + & + v_berr[12]^2 + & + v_berr[13]^2 + & + v_berr[14]^2 + & + v_berr[15]^2 + & + v_berr[16]^2 + & + v_berr[17]^2 + & + v_berr[18]^2 + & + v_berr[19]^2 + & + v_berr[20]^2 + & + v_berr[21]^2" + +variable bsq equal "c_born[1]^2 + & + c_born[2]^2 + & + c_born[3]^2 + & + c_born[4]^2 + & + c_born[5]^2 + & + c_born[6]^2 + & + c_born[7]^2 + & + c_born[8]^2 + & + c_born[9]^2 + & + c_born[10]^2 + & + c_born[11]^2 + & + c_born[12]^2 + & + c_born[13]^2 + & + c_born[14]^2 + & + c_born[15]^2 + & + c_born[16]^2 + & + c_born[17]^2 + & + c_born[18]^2 + & + c_born[19]^2 + & + c_born[20]^2 + & + c_born[21]^2" + +variable brelerr equal sqrt(v_berrsq/(v_bsq+${nugget})) + +thermo_style custom step temp pe etotal press v_frelerr v_vrelerr v_brelerr thermo ${nthermo} run ${nsteps} diff --git a/examples/numdiff/log.19Feb2022.log.numdiff.g++.1 b/examples/numdiff/log.19Feb2022.log.numdiff.g++.1 new file mode 100644 index 0000000000..34edf708bc --- /dev/null +++ b/examples/numdiff/log.19Feb2022.log.numdiff.g++.1 @@ -0,0 +1,197 @@ +LAMMPS (17 Feb 2022) +# Numerical difference calculation +# of error in forces, virial stress, and Born matrix + +# adjustable parameters + +variable nsteps index 500 # length of run +variable nthermo index 10 # thermo output interval +variable ndump index 500 # dump output interval +variable nlat index 3 # size of box +variable fdelta index 1.0e-4 # displacement size +variable vdelta index 1.0e-6 # strain size for numdiff/virial +variable bdelta index 1.0e-8 # strain size for numdiff Born matrix +variable temp index 10.0 # temperature +variable nugget equal 1.0e-6 # regularization for relerr + +units metal +atom_style atomic + +atom_modify map yes +lattice fcc 5.358000 +Lattice spacing in x,y,z = 5.358 5.358 5.358 +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 (16.074 16.074 16.074) + 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 (16.074 16.074 16.074) + create_atoms CPU = 0.000 seconds +mass 1 39.903 + +velocity all create ${temp} 2357 mom yes dist gaussian +velocity all create 10.0 2357 mom yes dist gaussian + +pair_style lj/cut 5.0 +pair_coeff 1 1 0.0102701 3.42 + +neighbor 0.0 bin +neigh_modify every 1 delay 0 check no + +timestep 0.001 +fix nve all nve + +# define numerical force calculation + +fix numforce all numdiff ${nthermo} ${fdelta} +fix numforce all numdiff 10 ${fdelta} +fix numforce all numdiff 10 1.0e-4 +variable ferrx atom f_numforce[1]-fx +variable ferry atom f_numforce[2]-fy +variable ferrz atom f_numforce[3]-fz +variable ferrsq atom v_ferrx^2+v_ferry^2+v_ferrz^2 +compute faverrsq all reduce ave v_ferrsq +variable fsq atom fx^2+fy^2+fz^2 +compute favsq all reduce ave v_fsq +variable frelerr equal sqrt(c_faverrsq/(c_favsq+${nugget})) +variable frelerr equal sqrt(c_faverrsq/(c_favsq+1e-06)) +dump errors all custom ${ndump} force_error.dump v_ferrx v_ferry v_ferrz +dump errors all custom 500 force_error.dump v_ferrx v_ferry v_ferrz + +# define numerical virial stress tensor calculation + +compute myvirial all pressure NULL virial +fix numvirial all numdiff/virial ${nthermo} ${vdelta} +fix numvirial all numdiff/virial 10 ${vdelta} +fix numvirial all numdiff/virial 10 1.0e-6 +variable errxx equal f_numvirial[1]-c_myvirial[1] +variable erryy equal f_numvirial[2]-c_myvirial[2] +variable errzz equal f_numvirial[3]-c_myvirial[3] +variable erryz equal f_numvirial[4]-c_myvirial[6] +variable errxz equal f_numvirial[5]-c_myvirial[5] +variable errxy equal f_numvirial[6]-c_myvirial[4] +variable verrsq equal "v_errxx^2 + v_erryy^2 + v_errzz^2 + v_erryz^2 + v_errxz^2 + v_errxy^2" +variable vsq equal "c_myvirial[1]^2 + c_myvirial[3]^2 + c_myvirial[3]^2 + c_myvirial[4]^2 + c_myvirial[5]^2 + c_myvirial[6]^2" +variable vrelerr equal sqrt(v_verrsq/(v_vsq+${nugget})) +variable vrelerr equal sqrt(v_verrsq/(v_vsq+1e-06)) + +# define numerical Born matrix calculation + +compute bornnum all born/matrix numdiff ${bdelta} myvirial +compute bornnum all born/matrix numdiff 1.0e-8 myvirial +compute born all born/matrix +variable berr vector c_bornnum-c_born +variable berrsq equal "v_berr[1]^2 + v_berr[2]^2 + v_berr[3]^2 + v_berr[4]^2 + v_berr[5]^2 + v_berr[6]^2 + v_berr[7]^2 + v_berr[8]^2 + v_berr[9]^2 + v_berr[10]^2 + v_berr[11]^2 + v_berr[12]^2 + v_berr[13]^2 + v_berr[14]^2 + v_berr[15]^2 + v_berr[16]^2 + v_berr[17]^2 + v_berr[18]^2 + v_berr[19]^2 + v_berr[20]^2 + v_berr[21]^2" + +variable bsq equal "c_born[1]^2 + c_born[2]^2 + c_born[3]^2 + c_born[4]^2 + c_born[5]^2 + c_born[6]^2 + c_born[7]^2 + c_born[8]^2 + c_born[9]^2 + c_born[10]^2 + c_born[11]^2 + c_born[12]^2 + c_born[13]^2 + c_born[14]^2 + c_born[15]^2 + c_born[16]^2 + c_born[17]^2 + c_born[18]^2 + c_born[19]^2 + c_born[20]^2 + c_born[21]^2" + +variable brelerr equal sqrt(v_berrsq/(v_bsq+${nugget})) +variable brelerr equal sqrt(v_berrsq/(v_bsq+1e-06)) + +thermo_style custom step temp pe etotal press v_frelerr v_vrelerr v_brelerr +thermo ${nthermo} +thermo 10 +run ${nsteps} +run 500 + 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 = 5 + ghost atom cutoff = 5 + binsize = 2.5, bins = 7 7 7 + 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) = 6.823 | 6.823 | 6.823 Mbytes +Step Temp PotEng TotEng Press v_frelerr v_vrelerr v_brelerr + 0 10 -6.6101864 -6.471878 947.70558 5.7012262e-09 1.4849734e-08 9.036398e-09 + 10 9.9357441 -6.6092976 -6.471878 949.3486 1.3828998e-08 1.9248385e-09 4.0233493e-09 + 20 9.7444561 -6.6066519 -6.4718779 954.23637 1.385204e-08 1.7476399e-09 4.0061081e-09 + 30 9.4311148 -6.6023181 -6.4718779 962.23331 1.4147226e-08 1.7647816e-09 3.3866543e-09 + 40 9.0043293 -6.5964152 -6.4718778 973.10762 1.4128155e-08 1.6390138e-09 3.2652821e-09 + 50 8.4762135 -6.5891108 -6.4718777 986.53572 1.4168048e-08 2.3910821e-09 4.7266627e-09 + 60 7.8621735 -6.5806179 -6.4718775 1002.1092 1.411958e-08 2.0683414e-09 2.6951001e-09 + 70 7.1805874 -6.5711908 -6.4718773 1019.3448 1.4139911e-08 1.6084571e-09 3.1477301e-09 + 80 6.4523557 -6.5611186 -6.4718771 1037.6974 1.4105096e-08 1.9929271e-09 3.4733802e-09 + 90 5.7003071 -6.5507169 -6.4718769 1056.5767 1.4084183e-08 1.750579e-09 4.310104e-09 + 100 4.9484503 -6.5403179 -6.4718767 1075.3674 1.4063796e-08 1.0250271e-09 2.9213594e-09 + 110 4.221081 -6.5302576 -6.4718765 1093.4526 1.400901e-08 1.389277e-09 4.3909721e-09 + 120 3.5417733 -6.520862 -6.4718763 1110.2388 1.4038158e-08 8.6231891e-10 2.5890696e-09 + 130 2.9323072 -6.5124324 -6.4718762 1125.183 1.4048645e-08 7.0840985e-10 3.388192e-09 + 140 2.411607 -6.5052306 -6.471876 1137.8182 1.3968429e-08 1.8508015e-09 3.2976031e-09 + 150 1.9947801 -6.4994654 -6.4718759 1147.7764 1.395965e-08 1.9484728e-09 4.2924605e-09 + 160 1.6923481 -6.4952825 -6.4718759 1154.8063 1.3948606e-08 1.5275137e-09 4.0204309e-09 + 170 1.5097515 -6.492757 -6.4718759 1158.7853 1.3845523e-08 1.5455e-09 4.8781309e-09 + 180 1.4471795 -6.4918916 -6.4718759 1159.7221 1.3788451e-08 1.578099e-09 3.0795316e-09 + 190 1.4997431 -6.4926187 -6.471876 1157.7529 1.374841e-08 2.142073e-09 2.4376961e-09 + 200 1.6579637 -6.4948072 -6.4718761 1153.1286 1.3674788e-08 2.111894e-09 3.7055708e-09 + 210 1.908522 -6.4982727 -6.4718763 1146.1965 1.3639408e-08 1.2386489e-09 3.160881e-09 + 220 2.23518 -6.5027908 -6.4718764 1137.3775 1.3524209e-08 1.7016573e-09 3.6982265e-09 + 230 2.6197892 -6.5081105 -6.4718766 1127.1415 1.3344007e-08 1.5843477e-09 3.7272821e-09 + 240 3.043298 -6.5139681 -6.4718768 1115.9815 1.3245227e-08 1.5502368e-09 3.898015e-09 + 250 3.4866901 -6.5201007 -6.4718769 1104.3906 1.3080142e-08 1.369987e-09 4.9133863e-09 + 260 3.9318061 -6.5262572 -6.4718771 1092.84 1.2885339e-08 1.0743728e-09 5.7271364e-09 + 270 4.3620216 -6.5322076 -6.4718772 1081.7617 1.2705966e-08 1.3618619e-09 2.3225062e-09 + 280 4.7627723 -6.5377504 -6.4718773 1071.5341 1.2480463e-08 1.4346869e-09 3.281167e-09 + 290 5.1219322 -6.542718 -6.4718774 1062.4716 1.2434727e-08 2.1935942e-09 2.8198924e-09 + 300 5.4300557 -6.5469796 -6.4718774 1054.8177 1.2321314e-08 8.2365886e-10 3.2731015e-09 + 310 5.6804997 -6.5504435 -6.4718774 1048.7409 1.2300884e-08 1.4855741e-09 4.1031988e-09 + 320 5.8694423 -6.5530567 -6.4718774 1044.3341 1.2483087e-08 1.8711589e-09 3.9368436e-09 + 330 5.9958115 -6.5548045 -6.4718774 1041.6165 1.2627617e-08 1.9256986e-09 4.3283764e-09 + 340 6.0611353 -6.555708 -6.4718774 1040.5369 1.2935701e-08 1.6609255e-09 3.8728039e-09 + 350 6.0693222 -6.5558211 -6.4718773 1040.9803 1.3218179e-08 1.985355e-09 2.618577e-09 + 360 6.0263776 -6.5552271 -6.4718773 1042.7755 1.3471701e-08 1.5125203e-09 2.936238e-09 + 370 5.9400629 -6.5540332 -6.4718772 1045.7049 1.3676495e-08 1.7364093e-09 2.9097362e-09 + 380 5.8195019 -6.5523657 -6.4718771 1049.515 1.3859995e-08 1.6834835e-09 2.7416302e-09 + 390 5.6747442 -6.5503635 -6.471877 1053.9288 1.3987553e-08 1.7893896e-09 2.8552537e-09 + 400 5.5162948 -6.5481719 -6.4718769 1058.6583 1.4091878e-08 1.4468098e-09 3.2733654e-09 + 410 5.3546269 -6.5459358 -6.4718768 1063.4182 1.4188438e-08 1.7231047e-09 3.3165187e-09 + 420 5.1996958 -6.5437929 -6.4718768 1067.9384 1.4205207e-08 1.3551982e-09 3.8687611e-09 + 430 5.0604771 -6.5418673 -6.4718767 1071.9767 1.4267199e-08 1.361845e-09 3.1210672e-09 + 440 4.9445529 -6.5402639 -6.4718766 1075.3292 1.4253464e-08 1.3945282e-09 2.6483572e-09 + 450 4.8577717 -6.5390637 -6.4718766 1077.8394 1.4240998e-08 1.8767323e-09 3.2040422e-09 + 460 4.8040023 -6.53832 -6.4718766 1079.4048 1.4242259e-08 1.4785379e-09 3.4402279e-09 + 470 4.7849977 -6.5380571 -6.4718766 1079.9795 1.4227939e-08 1.8623848e-09 4.3634918e-09 + 480 4.8003794 -6.5382699 -6.4718766 1079.5756 1.4215836e-08 1.2821795e-09 2.6846581e-09 + 490 4.8477405 -6.538925 -6.4718767 1078.2596 1.4186541e-08 2.47604e-09 3.2044632e-09 + 500 4.9228588 -6.539964 -6.4718767 1076.1469 1.4099819e-08 1.6653302e-09 3.267113e-09 +Loop time of 0.458483 on 1 procs for 500 steps with 108 atoms + +Performance: 94.224 ns/day, 0.255 hours/ns, 1090.552 timesteps/s +99.7% 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.0042278 | 0.0042278 | 0.0042278 | 0.0 | 0.92 +Neigh | 0.02481 | 0.02481 | 0.02481 | 0.0 | 5.41 +Comm | 0.002944 | 0.002944 | 0.002944 | 0.0 | 0.64 +Output | 0.014731 | 0.014731 | 0.014731 | 0.0 | 3.21 +Modify | 0.41122 | 0.41122 | 0.41122 | 0.0 | 89.69 +Other | | 0.0005545 | | | 0.12 + +Nlocal: 108 ave 108 max 108 min +Histogram: 1 0 0 0 0 0 0 0 0 0 +Nghost: 256 ave 256 max 256 min +Histogram: 1 0 0 0 0 0 0 0 0 0 +Neighs: 648 ave 648 max 648 min +Histogram: 1 0 0 0 0 0 0 0 0 0 + +Total # of neighbors = 648 +Ave neighs/atom = 6 +Neighbor list builds = 500 +Dangerous builds not checked +Total wall time: 0:00:00 diff --git a/examples/numdiff/log.19Feb2022.log.numdiff.g++.4 b/examples/numdiff/log.19Feb2022.log.numdiff.g++.4 new file mode 100644 index 0000000000..a692a6ad7f --- /dev/null +++ b/examples/numdiff/log.19Feb2022.log.numdiff.g++.4 @@ -0,0 +1,197 @@ +LAMMPS (17 Feb 2022) +# Numerical difference calculation +# of error in forces, virial stress, and Born matrix + +# adjustable parameters + +variable nsteps index 500 # length of run +variable nthermo index 10 # thermo output interval +variable ndump index 500 # dump output interval +variable nlat index 3 # size of box +variable fdelta index 1.0e-4 # displacement size +variable vdelta index 1.0e-6 # strain size for numdiff/virial +variable bdelta index 1.0e-8 # strain size for numdiff Born matrix +variable temp index 10.0 # temperature +variable nugget equal 1.0e-6 # regularization for relerr + +units metal +atom_style atomic + +atom_modify map yes +lattice fcc 5.358000 +Lattice spacing in x,y,z = 5.358 5.358 5.358 +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 (16.074 16.074 16.074) + 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 (16.074 16.074 16.074) + create_atoms CPU = 0.000 seconds +mass 1 39.903 + +velocity all create ${temp} 2357 mom yes dist gaussian +velocity all create 10.0 2357 mom yes dist gaussian + +pair_style lj/cut 5.0 +pair_coeff 1 1 0.0102701 3.42 + +neighbor 0.0 bin +neigh_modify every 1 delay 0 check no + +timestep 0.001 +fix nve all nve + +# define numerical force calculation + +fix numforce all numdiff ${nthermo} ${fdelta} +fix numforce all numdiff 10 ${fdelta} +fix numforce all numdiff 10 1.0e-4 +variable ferrx atom f_numforce[1]-fx +variable ferry atom f_numforce[2]-fy +variable ferrz atom f_numforce[3]-fz +variable ferrsq atom v_ferrx^2+v_ferry^2+v_ferrz^2 +compute faverrsq all reduce ave v_ferrsq +variable fsq atom fx^2+fy^2+fz^2 +compute favsq all reduce ave v_fsq +variable frelerr equal sqrt(c_faverrsq/(c_favsq+${nugget})) +variable frelerr equal sqrt(c_faverrsq/(c_favsq+1e-06)) +dump errors all custom ${ndump} force_error.dump v_ferrx v_ferry v_ferrz +dump errors all custom 500 force_error.dump v_ferrx v_ferry v_ferrz + +# define numerical virial stress tensor calculation + +compute myvirial all pressure NULL virial +fix numvirial all numdiff/virial ${nthermo} ${vdelta} +fix numvirial all numdiff/virial 10 ${vdelta} +fix numvirial all numdiff/virial 10 1.0e-6 +variable errxx equal f_numvirial[1]-c_myvirial[1] +variable erryy equal f_numvirial[2]-c_myvirial[2] +variable errzz equal f_numvirial[3]-c_myvirial[3] +variable erryz equal f_numvirial[4]-c_myvirial[6] +variable errxz equal f_numvirial[5]-c_myvirial[5] +variable errxy equal f_numvirial[6]-c_myvirial[4] +variable verrsq equal "v_errxx^2 + v_erryy^2 + v_errzz^2 + v_erryz^2 + v_errxz^2 + v_errxy^2" +variable vsq equal "c_myvirial[1]^2 + c_myvirial[3]^2 + c_myvirial[3]^2 + c_myvirial[4]^2 + c_myvirial[5]^2 + c_myvirial[6]^2" +variable vrelerr equal sqrt(v_verrsq/(v_vsq+${nugget})) +variable vrelerr equal sqrt(v_verrsq/(v_vsq+1e-06)) + +# define numerical Born matrix calculation + +compute bornnum all born/matrix numdiff ${bdelta} myvirial +compute bornnum all born/matrix numdiff 1.0e-8 myvirial +compute born all born/matrix +variable berr vector c_bornnum-c_born +variable berrsq equal "v_berr[1]^2 + v_berr[2]^2 + v_berr[3]^2 + v_berr[4]^2 + v_berr[5]^2 + v_berr[6]^2 + v_berr[7]^2 + v_berr[8]^2 + v_berr[9]^2 + v_berr[10]^2 + v_berr[11]^2 + v_berr[12]^2 + v_berr[13]^2 + v_berr[14]^2 + v_berr[15]^2 + v_berr[16]^2 + v_berr[17]^2 + v_berr[18]^2 + v_berr[19]^2 + v_berr[20]^2 + v_berr[21]^2" + +variable bsq equal "c_born[1]^2 + c_born[2]^2 + c_born[3]^2 + c_born[4]^2 + c_born[5]^2 + c_born[6]^2 + c_born[7]^2 + c_born[8]^2 + c_born[9]^2 + c_born[10]^2 + c_born[11]^2 + c_born[12]^2 + c_born[13]^2 + c_born[14]^2 + c_born[15]^2 + c_born[16]^2 + c_born[17]^2 + c_born[18]^2 + c_born[19]^2 + c_born[20]^2 + c_born[21]^2" + +variable brelerr equal sqrt(v_berrsq/(v_bsq+${nugget})) +variable brelerr equal sqrt(v_berrsq/(v_bsq+1e-06)) + +thermo_style custom step temp pe etotal press v_frelerr v_vrelerr v_brelerr +thermo ${nthermo} +thermo 10 +run ${nsteps} +run 500 + 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 = 5 + ghost atom cutoff = 5 + binsize = 2.5, bins = 7 7 7 + 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) = 6.816 | 6.816 | 6.816 Mbytes +Step Temp PotEng TotEng Press v_frelerr v_vrelerr v_brelerr + 0 10 -6.6101864 -6.471878 947.70558 1.9110624e-09 9.4407596e-10 3.1867416e-09 + 10 9.9369961 -6.6093149 -6.471878 949.31222 1.3055176e-08 4.996456e-10 2.7421655e-09 + 20 9.7500224 -6.6067289 -6.4718779 954.07898 1.3721178e-08 5.6039795e-10 2.3689718e-09 + 30 9.4448115 -6.6025075 -6.4718779 961.85502 1.3813156e-08 6.8451692e-10 1.9844663e-09 + 40 9.0305392 -6.5967776 -6.4718777 972.39819 1.3961749e-08 3.1134064e-10 1.7915052e-09 + 50 8.5196068 -6.5897109 -6.4718776 985.38158 1.3996941e-08 7.0149406e-10 2.002272e-09 + 60 7.9273388 -6.5815192 -6.4718775 1000.4024 1.4000005e-08 3.5766629e-10 2.4944703e-09 + 70 7.2715879 -6.5724494 -6.4718773 1016.9932 1.3996503e-08 6.2731503e-10 1.7010533e-09 + 80 6.5722375 -6.5627766 -6.4718771 1034.6361 1.3973603e-08 3.1142917e-10 2.808524e-09 + 90 5.8505991 -6.5527956 -6.4718769 1052.7794 1.3983301e-08 3.9931135e-10 2.6118214e-09 + 100 5.128708 -6.542811 -6.4718767 1070.8561 1.395586e-08 2.3152413e-10 2.8742755e-09 + 110 4.4285344 -6.5331269 -6.4718766 1088.305 1.3938374e-08 4.2173005e-10 2.3059886e-09 + 120 3.7711361 -6.5240343 -6.4718764 1104.5919 1.3915264e-08 2.5458038e-10 1.4864012e-09 + 130 3.1757964 -6.5158002 -6.4718762 1119.2319 1.3858843e-08 5.7490448e-10 2.6191823e-09 + 140 2.6591997 -6.5086551 -6.4718761 1131.8095 1.3814891e-08 3.5434633e-10 2.2009364e-09 + 150 2.2347034 -6.5027839 -6.471876 1141.9961 1.3781115e-08 5.0639594e-10 2.9032558e-09 + 160 1.9117661 -6.4983173 -6.471876 1149.564 1.3734288e-08 3.1954962e-10 2.6097446e-09 + 170 1.6955808 -6.4953273 -6.471876 1154.3946 1.3682252e-08 3.5426781e-10 2.9605676e-09 + 180 1.586949 -6.4938249 -6.471876 1156.4812 1.363e-08 4.0804881e-10 2.1707904e-09 + 190 1.5824056 -6.4937621 -6.4718761 1155.925 1.3532637e-08 4.0767685e-10 3.0091462e-09 + 200 1.6745831 -6.4950371 -6.4718762 1152.926 1.3455927e-08 2.953369e-10 2.5029298e-09 + 210 1.8527803 -6.4975018 -6.4718763 1147.7684 1.335224e-08 3.5042319e-10 3.0550064e-09 + 220 2.1036825 -6.5009721 -6.4718764 1140.8026 1.3239176e-08 3.5988448e-10 2.6852683e-09 + 230 2.4121721 -6.5052389 -6.4718766 1132.4243 1.3090019e-08 3.5004036e-10 2.8838602e-09 + 240 2.7621668 -6.5100798 -6.4718767 1123.0538 1.2946525e-08 4.1216361e-10 2.1105916e-09 + 250 3.1374274 -6.5152701 -6.4718768 1113.1152 1.277789e-08 5.9848318e-10 2.3087106e-09 + 260 3.5222906 -6.5205932 -6.471877 1103.0171 1.2591089e-08 2.0080182e-10 1.6969069e-09 + 270 3.9022942 -6.5258491 -6.4718771 1093.1369 1.2432232e-08 4.2494727e-10 1.7375594e-09 + 280 4.2646753 -6.5308612 -6.4718772 1083.8072 1.2268238e-08 6.1239266e-10 1.7005135e-09 + 290 4.598736 -6.5354816 -6.4718772 1075.306 1.2181179e-08 4.9338341e-10 2.1326848e-09 + 300 4.896078 -6.5395941 -6.4718773 1067.85 1.2098274e-08 3.4564838e-10 2.4199891e-09 + 310 5.150715 -6.543116 -6.4718773 1061.5918 1.2184958e-08 4.2383299e-10 2.2243759e-09 + 320 5.3590742 -6.5459978 -6.4718773 1056.6189 1.2312948e-08 3.5194185e-10 1.3856935e-09 + 330 5.5199009 -6.5482222 -6.4718773 1052.9565 1.2573918e-08 4.2401322e-10 2.9882e-09 + 340 5.6340787 -6.5498013 -6.4718773 1050.5719 1.2821551e-08 5.8802825e-10 2.7333289e-09 + 350 5.7043792 -6.5507736 -6.4718772 1049.3813 1.3067314e-08 4.0014945e-10 2.3564728e-09 + 360 5.7351548 -6.5511992 -6.4718772 1049.2581 1.331283e-08 4.1684815e-10 1.735621e-09 + 370 5.7319891 -6.5511553 -6.4718771 1050.042 1.354018e-08 3.8495426e-10 2.4460056e-09 + 380 5.7013193 -6.5507311 -6.4718771 1051.5496 1.3734888e-08 3.5333605e-10 2.5174342e-09 + 390 5.6500487 -6.5500219 -6.471877 1053.5847 1.3892287e-08 3.8154957e-10 1.77358e-09 + 400 5.5851679 -6.5491245 -6.471877 1055.9489 1.3988171e-08 5.8769536e-10 1.9262201e-09 + 410 5.5134009 -6.5481319 -6.4718769 1058.4508 1.4088779e-08 3.6754739e-10 2.7586362e-09 + 420 5.4408957 -6.547129 -6.4718769 1060.9152 1.4139924e-08 4.9030281e-10 3.2871245e-09 + 430 5.3729707 -6.5461895 -6.4718768 1063.1898 1.4173041e-08 5.2345074e-10 3.5995984e-09 + 440 5.3139284 -6.5453729 -6.4718768 1065.1506 1.4191516e-08 5.9481094e-10 2.5005297e-09 + 450 5.2669383 -6.5447229 -6.4718768 1066.7054 1.4168424e-08 3.0799668e-10 2.0864191e-09 + 460 5.2339881 -6.5442672 -6.4718768 1067.7958 1.4163444e-08 6.3927736e-10 2.2872669e-09 + 470 5.2158979 -6.544017 -6.4718768 1068.3968 1.413819e-08 5.5108262e-10 4.4334972e-09 + 480 5.2123873 -6.5439685 -6.4718768 1068.5155 1.4083227e-08 3.9249548e-10 2.5568235e-09 + 490 5.2221849 -6.544104 -6.4718768 1068.188 1.4035287e-08 2.1988631e-10 2.1264034e-09 + 500 5.2431716 -6.5443943 -6.4718768 1067.4759 1.3968666e-08 3.9100701e-10 3.290368e-09 +Loop time of 0.170182 on 4 procs for 500 steps with 108 atoms + +Performance: 253.846 ns/day, 0.095 hours/ns, 2938.035 timesteps/s +99.7% 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.0012069 | 0.0012994 | 0.0013512 | 0.2 | 0.76 +Neigh | 0.0048233 | 0.0050244 | 0.0053881 | 0.3 | 2.95 +Comm | 0.0072462 | 0.0078013 | 0.008175 | 0.4 | 4.58 +Output | 0.0080632 | 0.0081244 | 0.0082899 | 0.1 | 4.77 +Modify | 0.1476 | 0.14764 | 0.14768 | 0.0 | 86.75 +Other | | 0.0002961 | | | 0.17 + +Nlocal: 27 ave 31 max 24 min +Histogram: 1 0 1 0 1 0 0 0 0 1 +Nghost: 135 ave 138 max 131 min +Histogram: 1 0 0 0 0 1 0 1 0 1 +Neighs: 162 ave 191 max 148 min +Histogram: 1 2 0 0 0 0 0 0 0 1 + +Total # of neighbors = 648 +Ave neighs/atom = 6 +Neighbor list builds = 500 +Dangerous builds not checked +Total wall time: 0:00:00 diff --git a/examples/numdiff/log.28Jan2022.log.numdiff.g++.1 b/examples/numdiff/log.28Jan2022.log.numdiff.g++.1 deleted file mode 100644 index f32e72834f..0000000000 --- a/examples/numdiff/log.28Jan2022.log.numdiff.g++.1 +++ /dev/null @@ -1,175 +0,0 @@ -LAMMPS (7 Jan 2022) -# Numerical difference calculation -# of error in forces and virial stress - -# adjustable parameters - -variable nsteps index 500 # length of run -variable nthermo index 10 # thermo output interval -variable ndump index 500 # dump output interval -variable nlat index 3 # size of box -variable fdelta index 1.0e-4 # displacement size -variable vdelta index 1.0e-6 # strain size -variable temp index 10.0 # temperature - -units metal -atom_style atomic - -atom_modify map yes -lattice fcc 5.358000 -Lattice spacing in x,y,z = 5.358 5.358 5.358 -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 (16.074 16.074 16.074) - 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 (16.074 16.074 16.074) - create_atoms CPU = 0.000 seconds -mass 1 39.903 - -velocity all create ${temp} 2357 mom yes dist gaussian -velocity all create 10.0 2357 mom yes dist gaussian - -pair_style lj/cubic -pair_coeff * * 0.0102701 3.42 - -neighbor 0.0 bin -neigh_modify every 1 delay 0 check no - -timestep 0.001 -fix nve all nve - -# define numerical force calculation - -fix numforce all numdiff ${nthermo} ${fdelta} -fix numforce all numdiff 10 ${fdelta} -fix numforce all numdiff 10 1.0e-4 -variable ferrx atom f_numforce[1]-fx -variable ferry atom f_numforce[2]-fy -variable ferrz atom f_numforce[3]-fz -variable ferrsq atom v_ferrx^2+v_ferry^2+v_ferrz^2 -compute faverrsq all reduce ave v_ferrsq -variable fsq atom fx^2+fy^2+fz^2 -compute favsq all reduce ave v_fsq -variable frelerr equal sqrt(c_faverrsq/c_favsq) -dump errors all custom ${ndump} force_error.dump v_ferrx v_ferry v_ferrz -dump errors all custom 500 force_error.dump v_ferrx v_ferry v_ferrz - -# define numerical virial stress tensor calculation - -compute myvirial all pressure NULL virial -fix numvirial all numdiff/virial ${nthermo} ${vdelta} -fix numvirial all numdiff/virial 10 ${vdelta} -fix numvirial all numdiff/virial 10 1.0e-6 -variable errxx equal f_numvirial[1]-c_myvirial[1] -variable erryy equal f_numvirial[2]-c_myvirial[2] -variable errzz equal f_numvirial[3]-c_myvirial[3] -variable erryz equal f_numvirial[4]-c_myvirial[6] -variable errxz equal f_numvirial[5]-c_myvirial[5] -variable errxy equal f_numvirial[6]-c_myvirial[4] -variable verrsq equal "v_errxx^2 + v_erryy^2 + v_errzz^2 + v_erryz^2 + v_errxz^2 + v_errxy^2" -variable vsq equal "c_myvirial[1]^2 + c_myvirial[3]^2 + c_myvirial[3]^2 + c_myvirial[4]^2 + c_myvirial[5]^2 + c_myvirial[6]^2" -variable vrelerr equal sqrt(v_verrsq/v_vsq) - -thermo_style custom step temp pe etotal press v_frelerr v_vrelerr -thermo ${nthermo} -thermo 10 -run ${nsteps} -run 500 - 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 = 5.9407173 - ghost atom cutoff = 5.9407173 - binsize = 2.9703587, bins = 6 6 6 - 1 neighbor lists, perpetual/occasional/extra = 1 0 0 - (1) pair lj/cubic, perpetual - attributes: half, newton on - pair build: half/bin/atomonly/newton - stencil: half/bin/3d - bin: standard -Per MPI rank memory allocation (min/avg/max) = 6.083 | 6.083 | 6.083 Mbytes -Step Temp PotEng TotEng Press v_frelerr v_vrelerr - 0 10 -7.0259569 -6.8876486 28.564278 19203.344 1.5660292e-06 - 10 9.9376583 -7.0250947 -6.8876486 30.254762 1.5040965e-08 2.1991382e-07 - 20 9.7520139 -7.022527 -6.8876485 35.28505 1.4756358e-08 2.6265315e-06 - 30 9.4477557 -7.0183188 -6.8876485 43.519863 1.4688198e-08 2.6356166e-07 - 40 9.0330215 -7.0125826 -6.8876484 54.727797 1.4637921e-08 5.2292327e-08 - 50 8.5192918 -7.0054772 -6.8876483 68.585553 1.4587854e-08 7.1324716e-08 - 60 7.9212026 -6.997205 -6.8876481 84.684636 1.4525561e-08 3.1108149e-08 - 70 7.2562592 -6.9880081 -6.8876479 102.54088 1.450885e-08 3.2311094e-08 - 80 6.5444294 -6.9781627 -6.8876478 121.60715 1.4444738e-08 2.1776998e-08 - 90 5.8075961 -6.9679715 -6.8876476 141.2895 1.4493562e-08 2.0400898e-08 - 100 5.0688629 -6.957754 -6.8876474 160.9668 1.445455e-08 1.2636688e-08 - 110 4.3517145 -6.947835 -6.8876472 180.0135 1.4460371e-08 1.2528038e-08 - 120 3.6790589 -6.9385314 -6.887647 197.82486 1.4371757e-08 1.4489522e-08 - 130 3.0721984 -6.9301379 -6.8876468 213.84331 1.4364708e-08 1.2461922e-08 - 140 2.5497991 -6.9229125 -6.8876467 227.58429 1.4330926e-08 9.3913926e-09 - 150 2.1269443 -6.917064 -6.8876466 238.6596 1.4287002e-08 4.1510266e-09 - 160 1.8143642 -6.9127407 -6.8876465 246.79599 1.4282669e-08 7.7048281e-09 - 170 1.6179191 -6.9100237 -6.8876465 251.84748 1.42726e-08 1.2719973e-08 - 180 1.5383946 -6.9089239 -6.8876466 253.79991 1.4236534e-08 8.1200831e-09 - 190 1.5716287 -6.9093836 -6.8876467 252.76745 1.41706e-08 6.5670612e-09 - 200 1.7089493 -6.911283 -6.8876468 248.98142 1.4096463e-08 1.1685863e-08 - 210 1.9378716 -6.9144493 -6.8876469 242.77289 1.4008978e-08 1.1226902e-08 - 220 2.2429731 -6.9186692 -6.887647 234.55055 1.3886901e-08 9.9914102e-09 - 230 2.606862 -6.9237023 -6.8876472 224.77626 1.3864576e-08 1.1540228e-08 - 240 3.0111524 -6.9292941 -6.8876474 213.93996 1.3696314e-08 1.1697747e-08 - 250 3.4373794 -6.9351893 -6.8876475 202.53583 1.3626701e-08 1.0398197e-08 - 260 3.8678047 -6.9411426 -6.8876476 191.04084 1.3489489e-08 6.6603364e-09 - 270 4.2860853 -6.9469279 -6.8876478 179.89646 1.3312014e-08 1.1687917e-08 - 280 4.6777954 -6.9523457 -6.8876479 169.49404 1.3081144e-08 1.1336675e-08 - 290 5.030805 -6.9572282 -6.887648 160.16371 1.2947385e-08 1.7342825e-08 - 300 5.3355278 -6.9614428 -6.887648 152.16682 1.2893673e-08 1.7510534e-08 - 310 5.5850532 -6.964894 -6.887648 145.69148 1.2842022e-08 1.2782546e-08 - 320 5.7751794 -6.9675236 -6.8876481 140.85102 1.2903488e-08 1.5319437e-08 - 330 5.9043601 -6.9693103 -6.887648 137.68497 1.3076809e-08 1.1208999e-08 - 340 5.9735784 -6.9702676 -6.887648 136.16232 1.3296904e-08 1.891087e-08 - 350 5.9861549 -6.9704415 -6.887648 136.18679 1.3504051e-08 2.5783601e-08 - 360 5.947496 -6.9699067 -6.8876479 137.60397 1.3731112e-08 2.0556839e-08 - 370 5.8647874 -6.9687627 -6.8876478 140.2101 1.4009878e-08 2.1771736e-08 - 380 5.7466376 -6.9671285 -6.8876477 143.76234 1.4092054e-08 1.1085162e-08 - 390 5.6026773 -6.9651374 -6.8876477 147.99019 1.4282872e-08 2.0221602e-08 - 400 5.4431231 -6.9629305 -6.8876476 152.60787 1.4317739e-08 1.7076065e-08 - 410 5.2783192 -6.960651 -6.8876475 157.32722 1.4415075e-08 2.5031776e-08 - 420 5.1182723 -6.9584374 -6.8876474 161.87063 1.4441435e-08 2.2519289e-08 - 430 4.9722 -6.956417 -6.8876473 165.98344 1.4550624e-08 2.4512613e-08 - 440 4.8481153 -6.9547008 -6.8876473 169.44527 1.4544672e-08 1.4758301e-08 - 450 4.7524707 -6.9533779 -6.8876472 172.07964 1.4546492e-08 1.324687e-08 - 460 4.6898817 -6.9525122 -6.8876472 173.76132 1.4537475e-08 1.351367e-08 - 470 4.6629495 -6.9521397 -6.8876472 174.42109 1.4530458e-08 1.521106e-08 - 480 4.6721922 -6.9522675 -6.8876472 174.04742 1.4543785e-08 1.0905422e-08 - 490 4.7160887 -6.9528747 -6.8876473 172.68525 1.4545591e-08 2.0128525e-08 - 500 4.7912313 -6.953914 -6.8876473 170.43183 1.4438981e-08 1.6062775e-08 -Loop time of 0.837333 on 1 procs for 500 steps with 108 atoms - -Performance: 51.592 ns/day, 0.465 hours/ns, 597.134 timesteps/s -99.8% 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.0097726 | 0.0097726 | 0.0097726 | 0.0 | 1.17 -Neigh | 0.03095 | 0.03095 | 0.03095 | 0.0 | 3.70 -Comm | 0.005564 | 0.005564 | 0.005564 | 0.0 | 0.66 -Output | 0.0042451 | 0.0042451 | 0.0042451 | 0.0 | 0.51 -Modify | 0.78618 | 0.78618 | 0.78618 | 0.0 | 93.89 -Other | | 0.0006258 | | | 0.07 - -Nlocal: 108 ave 108 max 108 min -Histogram: 1 0 0 0 0 0 0 0 0 0 -Nghost: 558 ave 558 max 558 min -Histogram: 1 0 0 0 0 0 0 0 0 0 -Neighs: 972 ave 972 max 972 min -Histogram: 1 0 0 0 0 0 0 0 0 0 - -Total # of neighbors = 972 -Ave neighs/atom = 9 -Neighbor list builds = 500 -Dangerous builds not checked -Total wall time: 0:00:00 diff --git a/examples/numdiff/log.28Jan2022.log.numdiff.g++.4 b/examples/numdiff/log.28Jan2022.log.numdiff.g++.4 deleted file mode 100644 index fc56c4aee8..0000000000 --- a/examples/numdiff/log.28Jan2022.log.numdiff.g++.4 +++ /dev/null @@ -1,175 +0,0 @@ -LAMMPS (7 Jan 2022) -# Numerical difference calculation -# of error in forces and virial stress - -# adjustable parameters - -variable nsteps index 500 # length of run -variable nthermo index 10 # thermo output interval -variable ndump index 500 # dump output interval -variable nlat index 3 # size of box -variable fdelta index 1.0e-4 # displacement size -variable vdelta index 1.0e-6 # strain size -variable temp index 10.0 # temperature - -units metal -atom_style atomic - -atom_modify map yes -lattice fcc 5.358000 -Lattice spacing in x,y,z = 5.358 5.358 5.358 -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 (16.074 16.074 16.074) - 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 (16.074 16.074 16.074) - create_atoms CPU = 0.000 seconds -mass 1 39.903 - -velocity all create ${temp} 2357 mom yes dist gaussian -velocity all create 10.0 2357 mom yes dist gaussian - -pair_style lj/cubic -pair_coeff * * 0.0102701 3.42 - -neighbor 0.0 bin -neigh_modify every 1 delay 0 check no - -timestep 0.001 -fix nve all nve - -# define numerical force calculation - -fix numforce all numdiff ${nthermo} ${fdelta} -fix numforce all numdiff 10 ${fdelta} -fix numforce all numdiff 10 1.0e-4 -variable ferrx atom f_numforce[1]-fx -variable ferry atom f_numforce[2]-fy -variable ferrz atom f_numforce[3]-fz -variable ferrsq atom v_ferrx^2+v_ferry^2+v_ferrz^2 -compute faverrsq all reduce ave v_ferrsq -variable fsq atom fx^2+fy^2+fz^2 -compute favsq all reduce ave v_fsq -variable frelerr equal sqrt(c_faverrsq/c_favsq) -dump errors all custom ${ndump} force_error.dump v_ferrx v_ferry v_ferrz -dump errors all custom 500 force_error.dump v_ferrx v_ferry v_ferrz - -# define numerical virial stress tensor calculation - -compute myvirial all pressure NULL virial -fix numvirial all numdiff/virial ${nthermo} ${vdelta} -fix numvirial all numdiff/virial 10 ${vdelta} -fix numvirial all numdiff/virial 10 1.0e-6 -variable errxx equal f_numvirial[1]-c_myvirial[1] -variable erryy equal f_numvirial[2]-c_myvirial[2] -variable errzz equal f_numvirial[3]-c_myvirial[3] -variable erryz equal f_numvirial[4]-c_myvirial[6] -variable errxz equal f_numvirial[5]-c_myvirial[5] -variable errxy equal f_numvirial[6]-c_myvirial[4] -variable verrsq equal "v_errxx^2 + v_erryy^2 + v_errzz^2 + v_erryz^2 + v_errxz^2 + v_errxy^2" -variable vsq equal "c_myvirial[1]^2 + c_myvirial[3]^2 + c_myvirial[3]^2 + c_myvirial[4]^2 + c_myvirial[5]^2 + c_myvirial[6]^2" -variable vrelerr equal sqrt(v_verrsq/v_vsq) - -thermo_style custom step temp pe etotal press v_frelerr v_vrelerr -thermo ${nthermo} -thermo 10 -run ${nsteps} -run 500 - 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 = 5.9407173 - ghost atom cutoff = 5.9407173 - binsize = 2.9703587, bins = 6 6 6 - 1 neighbor lists, perpetual/occasional/extra = 1 0 0 - (1) pair lj/cubic, perpetual - attributes: half, newton on - pair build: half/bin/atomonly/newton - stencil: half/bin/3d - bin: standard -Per MPI rank memory allocation (min/avg/max) = 6.067 | 6.067 | 6.067 Mbytes -Step Temp PotEng TotEng Press v_frelerr v_vrelerr - 0 10 -7.0259569 -6.8876486 28.564278 10664.391 9.1481187e-08 - 10 9.9388179 -7.0251107 -6.8876486 30.21736 1.4771865e-08 1.3452512e-07 - 20 9.7572185 -7.022599 -6.8876485 35.123527 1.437525e-08 8.0966999e-07 - 30 9.4606673 -7.0184974 -6.8876484 43.132052 1.4375468e-08 1.990012e-08 - 40 9.0579092 -7.0129268 -6.8876483 54.000927 1.4350331e-08 1.7239028e-08 - 50 8.5607685 -7.0060508 -6.8876482 67.403151 1.4353284e-08 7.813181e-09 - 60 7.9838726 -6.9980717 -6.8876481 82.935358 1.4398078e-08 2.022316e-08 - 70 7.3442875 -6.9892255 -6.8876479 100.12892 1.434409e-08 7.5938179e-09 - 80 6.6610579 -6.9797757 -6.8876477 118.46358 1.4324787e-08 7.1972571e-09 - 90 5.9546462 -6.9700053 -6.8876476 137.38365 1.4322718e-08 4.3978378e-09 - 100 5.2462727 -6.9602077 -6.8876474 156.31651 1.4273172e-08 4.6728038e-09 - 110 4.5571706 -6.9506767 -6.8876472 174.69294 1.4266163e-08 3.522225e-09 - 120 3.9077807 -6.9416949 -6.887647 191.96859 1.42241e-08 3.5357511e-09 - 130 3.3169241 -6.9335227 -6.8876469 207.64566 1.4203813e-08 2.5182488e-09 - 140 2.8010028 -6.926387 -6.8876468 221.29333 1.4164215e-08 2.3112509e-09 - 150 2.3732854 -6.9204712 -6.8876467 232.5658 1.4134122e-08 1.9368963e-09 - 160 2.0433329 -6.9159076 -6.8876466 241.21647 1.4095473e-08 3.6806452e-09 - 170 1.8166184 -6.912772 -6.8876466 247.10715 1.4049531e-08 2.5319322e-09 - 180 1.6943727 -6.9110813 -6.8876467 250.21143 1.3997912e-08 1.952404e-09 - 190 1.6736731 -6.910795 -6.8876467 250.61203 1.3915487e-08 1.4271767e-09 - 200 1.7477659 -6.9118199 -6.8876468 248.49223 1.3850618e-08 2.4515718e-09 - 210 1.9065921 -6.9140167 -6.8876469 244.12226 1.3747916e-08 1.7957531e-09 - 220 2.1374676 -6.91721 -6.887647 237.84173 1.3674779e-08 2.6613511e-09 - 230 2.4258607 -6.9211989 -6.8876472 230.0395 1.3565712e-08 3.0777067e-09 - 240 2.7562034 -6.9257679 -6.8876473 221.13265 1.3440442e-08 1.7111501e-09 - 250 3.1126827 -6.9306984 -6.8876474 211.54566 1.3270914e-08 1.6690112e-09 - 260 3.4799641 -6.9357784 -6.8876476 201.69126 1.3105092e-08 2.1637558e-09 - 270 3.8438148 -6.9408108 -6.8876477 191.95361 1.2962846e-08 4.4613506e-09 - 280 4.191607 -6.9456212 -6.8876478 182.6745 1.2846383e-08 3.3730203e-09 - 290 4.5126922 -6.9500621 -6.8876478 174.14285 1.2710692e-08 2.2809889e-09 - 300 4.7986487 -6.9540172 -6.8876479 166.58747 1.2657778e-08 6.9880891e-09 - 310 5.0434083 -6.9574025 -6.8876479 160.17316 1.266381e-08 4.2486217e-09 - 320 5.243275 -6.9601668 -6.8876479 154.99974 1.279856e-08 5.1505673e-09 - 330 5.3968455 -6.9622908 -6.8876479 151.1038 1.2981831e-08 3.3339333e-09 - 340 5.5048468 -6.9637845 -6.8876479 148.46296 1.3159021e-08 1.7881393e-09 - 350 5.569902 -6.9646843 -6.8876479 147.00205 1.3439465e-08 5.6876219e-09 - 360 5.5962378 -6.9650485 -6.8876478 146.60113 1.3645943e-08 7.233847e-09 - 370 5.5893468 -6.9649531 -6.8876478 147.10471 1.3829581e-08 4.5514318e-09 - 380 5.5556199 -6.9644866 -6.8876477 148.33195 1.4005893e-08 4.2971846e-09 - 390 5.5019639 -6.9637444 -6.8876476 150.08725 1.4157454e-08 3.3564262e-09 - 400 5.4354239 -6.962824 -6.8876476 152.17073 1.4226422e-08 4.2393923e-09 - 410 5.3628267 -6.9618199 -6.8876475 154.38825 1.4302679e-08 3.8937698e-09 - 420 5.2904639 -6.960819 -6.8876475 156.56034 1.4381099e-08 4.315875e-09 - 430 5.2238282 -6.9598973 -6.8876474 158.52969 1.4426567e-08 4.2658185e-09 - 440 5.1674149 -6.9591171 -6.8876474 160.16704 1.4453381e-08 5.7055268e-09 - 450 5.1245913 -6.9585248 -6.8876474 161.37513 1.4449488e-08 4.4308801e-09 - 460 5.0975361 -6.9581506 -6.8876474 162.09077 1.4445596e-08 5.8269923e-09 - 470 5.0872416 -6.9580082 -6.8876474 162.28517 1.4444348e-08 4.8263194e-09 - 480 5.0935712 -6.9580957 -6.8876474 161.96268 1.4411666e-08 6.222228e-09 - 490 5.115362 -6.9583971 -6.8876474 161.15816 1.4369716e-08 3.3926077e-09 - 500 5.1505605 -6.958884 -6.8876474 159.9333 1.4288515e-08 3.8845251e-09 -Loop time of 0.252598 on 4 procs for 500 steps with 108 atoms - -Performance: 171.023 ns/day, 0.140 hours/ns, 1979.430 timesteps/s -99.8% 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.0021545 | 0.0022845 | 0.0023794 | 0.2 | 0.90 -Neigh | 0.0063887 | 0.0067604 | 0.0069752 | 0.3 | 2.68 -Comm | 0.01048 | 0.010916 | 0.011408 | 0.3 | 4.32 -Output | 0.0026603 | 0.0027399 | 0.0029738 | 0.3 | 1.08 -Modify | 0.2295 | 0.22952 | 0.22954 | 0.0 | 90.86 -Other | | 0.0003814 | | | 0.15 - -Nlocal: 27 ave 29 max 25 min -Histogram: 1 0 1 0 0 0 0 1 0 1 -Nghost: 325 ave 327 max 323 min -Histogram: 1 0 1 0 0 0 0 1 0 1 -Neighs: 243 ave 273 max 228 min -Histogram: 1 1 1 0 0 0 0 0 0 1 - -Total # of neighbors = 972 -Ave neighs/atom = 9 -Neighbor list builds = 500 -Dangerous builds not checked -Total wall time: 0:00:00 diff --git a/src/.gitignore b/src/.gitignore index 65077aa295..37e8b22893 100644 --- a/src/.gitignore +++ b/src/.gitignore @@ -453,6 +453,8 @@ /compute_basal_atom.h /compute_body_local.cpp /compute_body_local.h +/compute_born_matrix.cpp +/compute_born_matrix.h /compute_cnp_atom.cpp /compute_cnp_atom.h /compute_damage_atom.cpp diff --git a/src/EXTRA-COMPUTE/compute_born_matrix.cpp b/src/EXTRA-COMPUTE/compute_born_matrix.cpp new file mode 100644 index 0000000000..d0596e933f --- /dev/null +++ b/src/EXTRA-COMPUTE/compute_born_matrix.cpp @@ -0,0 +1,1221 @@ +/* ---------------------------------------------------------------------- + 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 "universe.h" +#include "update.h" + +#include +#include + +using namespace LAMMPS_NS; + +#define BIG 1000000000 +#define SMALL 1e-16 + +// this table is used to pick the 3d rij vector indices used to +// compute the 6 indices long Voigt stress vector + +static int constexpr sigma_albe[6][2] = { + {0, 0}, // s11 + {1, 1}, // s22 + {2, 2}, // s33 + {1, 2}, // s44 + {0, 2}, // s55 + {0, 1}, // s66 +}; + +// this table is used to pick the correct indices from the Voigt +// stress vector to compute the Cij matrix (21 terms, see doc) contribution + +static int constexpr C_albe[21][2] = { + {0, 0}, // C11 + {1, 1}, // C22 + {2, 2}, // C33 + {3, 3}, // C44 + {4, 4}, // C55 + {5, 5}, // 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 +}; + +// this table is used to pick the 3d rij vector indices used to +// compute the 21 indices long Cij matrix + +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), id_virial(nullptr), temp_x(nullptr), temp_f(nullptr) +{ + if (narg < 3) error->all(FLERR, "Illegal compute born/matrix command"); + + MPI_Comm_rank(world, &me); + + nvalues = 21; + + numflag = 0; + numdelta = 0.0; + + pairflag = bondflag = angleflag = dihedflag = impflag = 0; + if (narg == 3) { + pairflag = bondflag = angleflag = dihedflag = impflag = 1; + } else { + 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); + if (numdelta <= 0.0) error->all(FLERR, "Illegal compute born/matrix command"); + 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 if (strcmp(arg[iarg], "pair") == 0) { + pairflag = 1; + } else if (strcmp(arg[iarg], "bond") == 0) { + bondflag = 1; + } else if (strcmp(arg[iarg], "angle") == 0) { + angleflag = 1; + } else if (strcmp(arg[iarg], "dihedral") == 0) { + dihedflag = 1; + } else if (strcmp(arg[iarg], "improper") == 0) { + impflag = 1; + } else { + error->all(FLERR, "Illegal compute born/matrix command"); + } + ++iarg; + } + } + + if (pairflag) { + if (numflag) + error->all(FLERR, "Illegal compute born/matrix command: cannot mix numflag and other flags"); + if (force->pair) { + if (force->pair->born_matrix_enable == 0) + error->all(FLERR, "Pair style {} does not support compute born/matrix", force->pair_style); + } else { + pairflag = 0; + } + } + + if (bondflag) { + if (numflag) + error->all(FLERR, "Illegal compute born/matrix command: cannot mix numflag and other flags"); + if (force->bond) { + if (force->bond->born_matrix_enable == 0) + error->all(FLERR, "Bond style {} does not support compute born/matrix", force->bond_style); + } else { + bondflag = 0; + } + } + + if (angleflag) { + if (numflag) + error->all(FLERR, "Illegal compute born/matrix command: cannot mix numflag and other flags"); + if (force->angle) { + if (force->angle->born_matrix_enable == 0) + error->all(FLERR, "Angle style {} does not support compute born/matrix", + force->angle_style); + } else { + angleflag = 0; + } + } + + if (dihedflag) { + if (numflag) + error->all(FLERR, "Illegal compute born/matrix command: cannot mix numflag and other flags"); + if (force->dihedral) { + if (force->dihedral->born_matrix_enable == 0) + error->all(FLERR, "Dihedral style {} does not support compute born/matrix", + force->dihedral_style); + } else { + dihedflag = 0; + } + } + + if (impflag) { + if (numflag) + error->all(FLERR, "Illegal compute born/matrix command: cannot mix numflag and other flags"); + if (force->improper) { + if (force->improper->born_matrix_enable == 0) + error->all(FLERR, "Improper style {} does not support compute born/matrix", + force->improper_style); + } else { + impflag = 0; + } + } + + if (force->kspace) { + if (!numflag && (comm->me == 0)) + error->all(FLERR, "KSpace contribution not supported by compute born/matrix"); + } + + // 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); + delete[] id_virial; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputeBornMatrix::init() +{ + dt = update->dt; + + if (!numflag) { + + // need an occasional half neighbor list + + neighbor->add_request(this, NeighConst::REQ_FULL | NeighConst::REQ_OCCASIONAL); + + } 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 + // This table is used for consistency between numdiff and analytical + // ordering of the terms. + + for (int m = 0; m < nvalues; m++) { + int a = C_albe[m][0]; + int b = C_albe[m][1]; + revalbe[a][b] = m; + revalbe[b][a] = m; + } + + // reorder LAMMPS virial vector to Voigt order + + virialVtoV[0] = 0; + virialVtoV[1] = 1; + virialVtoV[2] = 2; + virialVtoV[3] = 5; + virialVtoV[4] = 4; + virialVtoV[5] = 3; + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputeBornMatrix::init_list(int /* id */, NeighList *ptr) +{ + list = ptr; +} + +/* ---------------------------------------------------------------------- + compute output vector +------------------------------------------------------------------------- */ + +void ComputeBornMatrix::compute_vector() +{ + invoked_vector = 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(); + 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(); + + // convert from pressure to energy units + + double inv_nktv2p = 1.0 / force->nktv2p; + double volume = domain->xprd * domain->yprd * domain->zprd; + for (int m = 0; m < nvalues; m++) { values_global[m] *= inv_nktv2p * volume; } + } + + 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; + + // Add contribution to Born tensor + + pair_pref = dupair = du2pair = 0.0; + pair->born_matrix(i, j, itype, jtype, rsq, factor_coul, factor_lj, dupair, du2pair); + pair_pref = 0.5 * du2pair - dupair * rinv; + + // See albemunu in compute_born_matrix.h for indices order. + + a = b = c = 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 vec_index; + + // grow arrays if necessary + + int nall = atom->nlocal + atom->nghost; + if (nall > maxatom) reallocate(); + + // store copy of current forces for owned and ghost atoms + + double **x = atom->x; + double **f = atom->f; + + 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 stress finite difference in each direction + + int flag, allflag; + + for (int idir = 0; idir < NDIR_VIRIAL; idir++) { + + // forward + + displace_atoms(nall, idir, 1.0); + force_clear(nall); + update_virial(); + for (int jdir = 0; jdir < NDIR_VIRIAL; jdir++) { + vec_index = revalbe[idir][jdir]; + values_global[vec_index] = compute_virial->vector[virialVtoV[jdir]]; + } + restore_atoms(nall, idir); + + // backward + + displace_atoms(nall, idir, -1.0); + force_clear(nall); + update_virial(); + for (int jdir = 0; jdir < NDIR_VIRIAL; jdir++) { + vec_index = revalbe[idir][jdir]; + values_global[vec_index] -= compute_virial->vector[virialVtoV[jdir]]; + } + restore_atoms(nall, idir); + } + + // apply derivative factor + + double denominator = -0.5 / numdelta; + for (int m = 0; m < nvalues; m++) values_global[m] *= denominator; + + // recompute virial so all virial and energy contributions are as before + // also needed for virial stress addon contributions to Born matrix + + update_virial(); + + // add on virial terms + + virial_addon(); + + // 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]; +} + +/* ---------------------------------------------------------------------- + displace position of all owned and ghost atoms +------------------------------------------------------------------------- */ + +void ComputeBornMatrix::displace_atoms(int nall, int idir, double magnitude) +{ + double **x = atom->x; + + // NOTE: virial_addon() expressions predicated on + // shear strain fields (l != k) being symmetric here + int k = dirlist[idir][0]; + int l = dirlist[idir][1]; + + // axial strain + + if (l == k) + for (int i = 0; i < nall; i++) + x[i][k] = temp_x[i][k] + numdelta * magnitude * (temp_x[i][l] - fixedpoint[l]); + + // symmetric shear strain + + else + for (int i = 0; i < nall; i++) { + x[i][k] = temp_x[i][k] + 0.5 * numdelta * magnitude * (temp_x[i][l] - fixedpoint[l]); + x[i][l] = temp_x[i][l] + 0.5 * numdelta * magnitude * (temp_x[i][k] - fixedpoint[k]); + } +} + +/* ---------------------------------------------------------------------- + restore position of all owned and ghost atoms +------------------------------------------------------------------------- */ + +void ComputeBornMatrix::restore_atoms(int nall, int idir) +{ + + // reset only idir coord + + int k = dirlist[idir][0]; + int l = dirlist[idir][1]; + double **x = atom->x; + if (l == k) + for (int i = 0; i < nall; i++) x[i][k] = temp_x[i][k]; + else + for (int i = 0; i < nall; i++) { + x[i][l] = temp_x[i][l]; + x[i][k] = temp_x[i][k]; + } +} + +/* ---------------------------------------------------------------------- + evaluate potential forces and virial + same logic as in Verlet +------------------------------------------------------------------------- */ + +void ComputeBornMatrix::update_virial() +{ + int eflag = 0; + + // this may not be completely general + // but it works for lj/cut and sw pair styles + // and compute stress/atom output is unaffected + + int vflag = VIRIAL_PAIR; + + if (force->pair) 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 (force->kspace) force->kspace->compute(eflag, vflag); + + compute_virial->compute_vector(); +} + +/* ---------------------------------------------------------------------- + calculate virial stress addon terms to the Born matrix + this is based on original code of Dr. Yubao Zhen + described here: Comp. Phys. Comm. 183 (2012) 261-265 + as well as Yoshimoto et al., PRB, 71 (2005) 184108, Eq 15.and eq A3. +------------------------------------------------------------------------- */ + +void ComputeBornMatrix::virial_addon() +{ + + int kd, nd, id, jd; + int m; + + double *sigv = compute_virial->vector; + + // This way of doing is not very elegant but is correct. + // The complete Cijkl terms are the sum of symmetric terms + // computed in compute_numdiff and virial stress terms. + // The viral terms are not symmetric in the tensor computation. + // For example: + // C2323 = s33+D2323; C3232 = s22+D3232 etc... + // However there are two symmetry breaking when reducing + // the 4-rank tensor to a 2-rank tensor + // Cijkl = (Bijkl+Bjikl+Bijlk+Bjilk)/4. = (Bijkl+Bjilk)/2. + // and when computing only the 21 independant term. + + // these expressions have been verified + // correct to about 1e-7 compared + // to the analytic expressions for lj/cut, + // predicated on shear strain fields being + // symmetric in displace_atoms() + + values_global[0] += 2.0 * sigv[0]; + values_global[1] += 2.0 * sigv[1]; + values_global[2] += 2.0 * sigv[2]; + + values_global[3] += 0.5 * (sigv[1] + sigv[2]); + values_global[4] += 0.5 * (sigv[0] + sigv[2]); + values_global[5] += 0.5 * (sigv[0] + sigv[1]); + values_global[6] += 0.0; + values_global[7] += 0.0; + values_global[8] += 0.0; + values_global[9] += sigv[4]; + values_global[10] += sigv[3]; + values_global[11] += 0.0; + values_global[12] += sigv[5]; + values_global[13] += 0.0; + values_global[14] += sigv[3]; + values_global[15] += sigv[5]; + values_global[16] += sigv[4]; + values_global[17] += 0.0; + values_global[18] += 0.5 * sigv[3]; + values_global[19] += 0.5 * sigv[4]; + values_global[20] += 0.5 * sigv[5]; +} + +/* ---------------------------------------------------------------------- + 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 +---------------------------------------------------------------------- */ + +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 + + for (atom1 = 0; atom1 < nlocal; atom1++) { + if (!(mask[atom1] & groupbit)) continue; + + if (molecular == 1) + nb = num_bond[atom1]; + else { + if (molindex[atom1] < 0) continue; + imol = molindex[atom1]; + iatom = molatom[atom1]; + nb = onemols[imol]->num_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 = dupair = du2pair = 0.0; + bond->born_matrix(btype, rsq, atom1, atom2, dupair, du2pair); + pair_pref = du2pair - dupair * rinv; + + a = b = c = d = 0; + for (m = 0; m < 21; 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; + } + } + } +} + +/* ---------------------------------------------------------------------- + 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 +---------------------------------------------------------------------- */ + +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, 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; + double d2lncos; + + // Initializing array for intermediate cos derivatives + // w regard to strain + for (i = 0; i < 6; i++) dcos[i] = 0; + + 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); + + 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); + + r1r2 = delx1 * delx2 + dely1 * dely2 + delz1 * delz2; + r1r2inv = 0; + if (r1r2 == 0.) { + // Worst case scenario. A 0 cosine has an undefined logarithm. + // We then add a small amount of the third bond to the first one + // so they are not orthogonal anymore and recompute. + del1[0] += SMALL * del2[0]; + del1[1] += SMALL * del2[1]; + del1[2] += SMALL * del2[2]; + r1r2 = del1[0] * del2[0] + del1[1] * del2[1] + del1[2] * del2[2]; + r1r2inv = 1 / r1r2; + + // cost = cosine of angle + cost = r1r2 / (r1 * r2); + + } else { + r1r2inv = 1 / r1r2; + cost = r1r2 / (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 = b = c = d = 0; + // clang-format off + for (m = 0; m<6; m++) { + a = sigma_albe[m][0]; + b = sigma_albe[m][1]; + dcos[m] = cost*((del1[a]*del2[b]+del1[b]*del2[a])*r1r2inv - + del1[a]*del1[b]*rsq1inv - del2[a]*del2[b]*rsq2inv); + } + for (m = 0; m<21; m++) { + a = albemunu[m][0]; + b = albemunu[m][1]; + c = albemunu[m][2]; + d = albemunu[m][3]; + e = C_albe[m][0]; + f = C_albe[m][1]; + d2lncos = 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 = cost*d2lncos + dcos[e]*dcos[f]*cinv; + values_local[m] += duang*d2cos + du2ang*dcos[e]*dcos[f]; + } + // clang-format on + } + } +} + +/* ---------------------------------------------------------------------- + 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 +------------------------------------------------------------------------- */ + +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, dotab, rgsq, rg, ra2inv, rb2inv, dotabinv, 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 al, be, mu, nu, e, f; + double b1sq; + double b2sq; + double b3sq; + double b1b2; + double b1b3; + double b2b3; + double b1[3]; + double b2[3]; + double b3[3]; + + // 1st and 2nd order derivatives of the dot products + double dab[6]; + double daa[6]; + double dbb[6]; + double d2ab; + double d2aa; + double d2bb; + + double dcos[6]; + double d2cos; + + for (i = 0; i < 6; i++) dab[i] = daa[i] = dbb[i] = dcos[i] = 0; + + for (atom2 = 0; atom2 < nlocal; atom2++) { + if (!(mask[atom2] & groupbit)) continue; + + // if (molecular == 1) + // nd = num_dihedral[atom2]; + if (molecular == Atom::MOLECULAR) + 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(i, atom1, atom2, atom3, atom4, dudih, du2dih); + + vb1x = x[atom2][0] - x[atom1][0]; + vb1y = x[atom2][1] - x[atom1][1]; + vb1z = x[atom2][2] - x[atom1][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]; + + 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]; + + // a = b1xb2; b = b2xb3 + ax = vb1y * vb2z - vb1z * vb2y; + ay = vb1z * vb2x - vb1x * vb2z; + az = vb1x * vb2y - vb1y * vb2x; + bx = vb2y * vb3z - vb2z * vb3y; + by = vb2z * vb3x - vb2x * vb3z; + bz = vb2x * vb3y - vb2y * vb3x; + + rasq = ax * ax + ay * ay + az * az; + rbsq = bx * bx + by * by + bz * bz; + rgsq = vb2x * vb2x + vb2y * vb2y + vb2z * vb2z; + dotab = ax * bx + ay * by + az * bz; + rg = sqrt(rgsq); + + ra2inv = rb2inv = rabinv = dotabinv = 0.0; + if (rasq > 0) ra2inv = 1.0 / rasq; + if (rbsq > 0) rb2inv = 1.0 / rbsq; + if (dotab != 0.) dotabinv = 1.0 / dotab; + rabinv = sqrt(ra2inv * rb2inv); + + co = (ax * bx + ay * by + az * bz) * rabinv; + si = sqrt(b2sq) * rabinv * (ax * vb3x + ay * vb3y + az * vb3z); + + if (co == 0.) { + // Worst case scenario. A 0 cosine has an undefined logarithm. + // We then add a small amount of the third bond to the first one + // so they are not orthogonal anymore and recompute. + b1[0] += SMALL * b3[0]; + b1[1] += SMALL * b3[1]; + b1[2] += SMALL * b3[2]; + b1sq = b1[0] * b1[0] + b1[1] * b1[1] + b1[2] * b1[2]; + 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 = b1[1] * b2[2] - b1[2] * b2[1]; + ay = b1[2] * b2[0] - b1[0] * b2[2]; + az = b1[0] * b2[1] - b1[1] * b2[0]; + bx = b2[1] * b3[2] - b2[2] * b3[1]; + by = b2[2] * b3[0] - b2[0] * b3[2]; + bz = b2[0] * b3[1] - b2[1] * b3[0]; + rasq = ax * ax + ay * ay + az * az; + rbsq = bx * bx + by * by + bz * bz; + dotab = ax * bx + ay * by + az * bz; + ra2inv = rb2inv = rabinv = dotabinv = 0.0; + if (rasq > 0) ra2inv = 1.0 / rasq; + if (rbsq > 0) rb2inv = 1.0 / rbsq; + if (dotab != 0.) dotabinv = 1.0 / dotab; + rabinv = sqrt(ra2inv * rb2inv); + co = (ax * bx + ay * by + az * bz) * rabinv; + si = sqrt(b2sq) * 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); + + al = be = mu = nu = e = f = 0; + // clang-format off + for (m = 0; m<6; m++) { + al = sigma_albe[m][0]; + be = sigma_albe[m][1]; + + daa[m] = 2*(b2sq*b1[al]*b1[be] + + b1sq*b2[al]*b2[be] + - b1b2*(b1[al]*b2[be]+b2[al]*b1[be])); + + dbb[m] = 2*(b2sq*b3[al]*b3[be] + + b3sq*b2[al]*b2[be] + - b2b3*(b3[al]*b2[be]+b2[al]*b3[be])); + + dab[m] = b1b2*(b2[al]*b3[be]+b3[al]*b2[be]) + + b2b3*(b1[al]*b2[be]+b2[al]*b1[be]) + - b1b3*(b2[al]*b2[be]+b2[al]*b2[be]) + - b2sq*(b1[al]*b3[be]+b3[al]*b1[be]); + + + dcos[m] = 0.5*co*(2*dotabinv*dab[m] - ra2inv*daa[m] - rb2inv*dbb[m]); + } + + for (m = 0; m<21; m++) { + al = albemunu[m][0]; + be = albemunu[m][1]; + mu = albemunu[m][2]; + nu = albemunu[m][3]; + e = C_albe[m][0]; + f = C_albe[m][1]; + + d2aa = 4*(b1[al]*b1[be]*b2[mu]*b2[nu] + b1[mu]*b1[nu]*b2[al]*b2[be]) + - 2*(b1[al]*b2[be]+b1[be]*b2[al])*(b1[mu]*b2[nu]+b1[nu]*b2[mu]); + + d2bb = 4*(b3[al]*b3[be]*b2[mu]*b2[nu] + b3[mu]*b3[nu]*b2[al]*b2[be]) + - 2*(b3[al]*b2[be]+b3[be]*b2[al])*(b3[mu]*b2[nu]+b3[nu]*b2[mu]); + + d2ab = (b1[al]*b2[be]+b2[al]*b1[be])*(b3[mu]*b2[nu]+b2[mu]*b3[nu]) + + (b3[al]*b2[be]+b2[al]*b3[be])*(b1[mu]*b2[nu]+b2[mu]*b1[nu]) + - (b1[al]*b3[be]+b3[al]*b1[be])*(b2[mu]*b2[nu]+b2[mu]*b2[nu]) + - (b2[al]*b2[be]+b2[al]*b2[be])*(b1[mu]*b3[nu]+b3[mu]*b1[nu]); + + d2cos = -2*dotabinv*dotabinv*dab[e]*dab[f] + 2*dotabinv*d2ab + + ra2inv*ra2inv*daa[e]*daa[f] - ra2inv*d2aa + + rb2inv*rb2inv*dbb[e]*dbb[f] - rb2inv*d2bb + + 2*dcos[e]*dcos[f]/co/co; + d2cos *= 0.5*co; + + values_local[m] += dudih*d2cos; + values_local[m] += du2dih*dcos[e]*dcos[f]; + } + // clang-format on + } + } +} diff --git a/src/EXTRA-COMPUTE/compute_born_matrix.h b/src/EXTRA-COMPUTE/compute_born_matrix.h new file mode 100644 index 0000000000..9991ebcc69 --- /dev/null +++ b/src/EXTRA-COMPUTE/compute_born_matrix.h @@ -0,0 +1,97 @@ +/* ---------------------------------------------------------------------- + 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), Aidan Thompson (Sandia) + --------------------------------------------------------------------------*/ + +#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 **); + ~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 virialVtoV[NDIR_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 +}; +} // namespace LAMMPS_NS + +#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/EXTRA-MOLECULE/dihedral_nharmonic.cpp b/src/EXTRA-MOLECULE/dihedral_nharmonic.cpp index 41a3d66ba6..fc1717bc94 100644 --- a/src/EXTRA-MOLECULE/dihedral_nharmonic.cpp +++ b/src/EXTRA-MOLECULE/dihedral_nharmonic.cpp @@ -25,6 +25,7 @@ #include "force.h" #include "memory.h" #include "neighbor.h" +#include "domain.h" #include @@ -39,6 +40,7 @@ DihedralNHarmonic::DihedralNHarmonic(LAMMPS *lmp) : Dihedral(lmp) { writedata = 1; a = nullptr; + born_matrix_enable = 1; } /* ---------------------------------------------------------------------- */ @@ -336,3 +338,64 @@ void DihedralNHarmonic::write_data(FILE *fp) } } + +/* ----------------------------------------------------------------------*/ + +void DihedralNHarmonic::born_matrix(int nd, int i1, int i2, int i3, int i4, + double &dudih, double &du2dih) { + int i,type; + double vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z,vb2xm,vb2ym,vb2zm; + double ax,ay,az,bx,by,bz,rasq,rbsq,rgsq,rg,rginv,ra2inv,rb2inv,rabinv; + double c,s,kf; + + int **dihedrallist = neighbor->dihedrallist; + double **x = atom->x; + + int ndihedrallist = neighbor->ndihedrallist; + type = dihedrallist[nd][4]; + + vb1x = x[i1][0] - x[i2][0]; + vb1y = x[i1][1] - x[i2][1]; + vb1z = x[i1][2] - x[i2][2]; + + vb2x = x[i3][0] - x[i2][0]; + vb2y = x[i3][1] - x[i2][1]; + vb2z = x[i3][2] - x[i2][2]; + + vb2xm = -vb2x; + vb2ym = -vb2y; + vb2zm = -vb2z; + + vb3x = x[i4][0] - x[i3][0]; + vb3y = x[i4][1] - x[i3][1]; + vb3z = x[i4][2] - x[i3][2]; + + // c,s calculation + + 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; + + ra2inv = rb2inv = 0.0; + if (rasq > 0) ra2inv = 1.0/rasq; + if (rbsq > 0) rb2inv = 1.0/rbsq; + rabinv = sqrt(ra2inv*rb2inv); + + c = (ax*bx + ay*by + az*bz)*rabinv; + + dudih = 0.0; + du2dih = 0.0; + for (i = 1; i < nterms[type]; i++) { + dudih += this->a[type][i]*i*pow(c,i-1); + } + + for (i = 2; i < nterms[type]; i++) { + du2dih += this->a[type][i]*i*(i-1)*pow(c, i-2); + } +} diff --git a/src/EXTRA-MOLECULE/dihedral_nharmonic.h b/src/EXTRA-MOLECULE/dihedral_nharmonic.h index 184cf1f740..28dd2e084e 100644 --- a/src/EXTRA-MOLECULE/dihedral_nharmonic.h +++ b/src/EXTRA-MOLECULE/dihedral_nharmonic.h @@ -33,6 +33,7 @@ class DihedralNHarmonic : public Dihedral { void write_restart(FILE *) override; void read_restart(FILE *) override; void write_data(FILE *) override; + void born_matrix(int /*dtype*/, int, int, int, int, double &, double &) override; protected: int *nterms; diff --git a/src/MOLECULE/angle_cosine.cpp b/src/MOLECULE/angle_cosine.cpp index 260ebc3948..f11fff8af0 100644 --- a/src/MOLECULE/angle_cosine.cpp +++ b/src/MOLECULE/angle_cosine.cpp @@ -29,7 +29,10 @@ using MathConst::MY_PI; /* ---------------------------------------------------------------------- */ -AngleCosine::AngleCosine(LAMMPS *_lmp) : Angle(_lmp) {} +AngleCosine::AngleCosine(LAMMPS *_lmp) : Angle(_lmp) +{ + born_matrix_enable = 1; +} /* ---------------------------------------------------------------------- */ @@ -235,6 +238,14 @@ double AngleCosine::single(int type, int i1, int i2, int i3) return k[type] * (1.0 + c); } +/* ---------------------------------------------------------------------- */ + +void AngleCosine::born_matrix(int type, int i1, int i2, int i3, double &du, double &du2) +{ + du2 = 0; + du = k[type]; +} + /* ---------------------------------------------------------------------- return ptr to internal members upon request ------------------------------------------------------------------------ */ diff --git a/src/MOLECULE/angle_cosine.h b/src/MOLECULE/angle_cosine.h index e249e7a44c..d86e0dc92f 100644 --- a/src/MOLECULE/angle_cosine.h +++ b/src/MOLECULE/angle_cosine.h @@ -35,6 +35,7 @@ class AngleCosine : public Angle { void read_restart(FILE *) override; void write_data(FILE *) override; double single(int, int, int, int) override; + void born_matrix(int type, int i1, int i2, int i3, double& du, double& du2) override; void *extract(const char *, int &) override; protected: diff --git a/src/MOLECULE/angle_cosine_squared.cpp b/src/MOLECULE/angle_cosine_squared.cpp index 262f345118..32e6f39b36 100644 --- a/src/MOLECULE/angle_cosine_squared.cpp +++ b/src/MOLECULE/angle_cosine_squared.cpp @@ -38,6 +38,7 @@ AngleCosineSquared::AngleCosineSquared(LAMMPS *_lmp) : Angle(_lmp) { k = nullptr; theta0 = nullptr; + born_matrix_enable = 1; } /* ---------------------------------------------------------------------- */ @@ -262,3 +263,32 @@ double AngleCosineSquared::single(int type, int i1, int i2, int i3) double tk = k[type] * dcostheta; return tk * dcostheta; } + +/* ---------------------------------------------------------------------- */ + +void AngleCosineSquared::born_matrix(int type, int i1, int i2, int i3, double& du, double& du2) +{ + double **x = atom->x; + + double delx1 = x[i1][0] - x[i2][0]; + double dely1 = x[i1][1] - x[i2][1]; + double delz1 = x[i1][2] - x[i2][2]; + domain->minimum_image(delx1,dely1,delz1); + double r1 = sqrt(delx1*delx1 + dely1*dely1 + delz1*delz1); + + double delx2 = x[i3][0] - x[i2][0]; + double dely2 = x[i3][1] - x[i2][1]; + double delz2 = x[i3][2] - x[i2][2]; + domain->minimum_image(delx2,dely2,delz2); + double r2 = sqrt(delx2*delx2 + dely2*dely2 + delz2*delz2); + + double c = delx1*delx2 + dely1*dely2 + delz1*delz2; + c /= r1*r2; + if (c > 1.0) c = 1.0; + if (c < -1.0) c = -1.0; + + double dcostheta = c - cos(theta0[type]); + double tk = k[type] * dcostheta; + du2 = 2*k[type]; + du = du2*dcostheta; +} diff --git a/src/MOLECULE/angle_cosine_squared.h b/src/MOLECULE/angle_cosine_squared.h index d8f282ac47..f2800c05b1 100644 --- a/src/MOLECULE/angle_cosine_squared.h +++ b/src/MOLECULE/angle_cosine_squared.h @@ -35,6 +35,7 @@ class AngleCosineSquared : public Angle { void read_restart(FILE *) override; void write_data(FILE *) override; double single(int, int, int, int) override; + void born_matrix(int type, int i1, int i2, int i3, double& du, double& du2) override; protected: double *k, *theta0; diff --git a/src/MOLECULE/bond_harmonic.cpp b/src/MOLECULE/bond_harmonic.cpp index 687a871f17..c946d9ff89 100644 --- a/src/MOLECULE/bond_harmonic.cpp +++ b/src/MOLECULE/bond_harmonic.cpp @@ -27,7 +27,10 @@ using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ -BondHarmonic::BondHarmonic(LAMMPS *_lmp) : Bond(_lmp) {} +BondHarmonic::BondHarmonic(LAMMPS *_lmp) : Bond(_lmp) +{ + born_matrix_enable = 1; +} /* ---------------------------------------------------------------------- */ @@ -197,6 +200,20 @@ double BondHarmonic::single(int type, double rsq, int /*i*/, int /*j*/, double & return rk * dr; } + +/* ---------------------------------------------------------------------- */ + +void BondHarmonic::born_matrix(int type, double rsq, int /*i*/, int /*j*/, + double &du, double& du2) +{ + double r = sqrt(rsq); + double dr = r - r0[type]; + du2 = 0.0; + du = 0.0; + du2 = 2*k[type]; + if (r > 0.0) du = du2*dr; +} + /* ---------------------------------------------------------------------- return ptr to internal members upon request ------------------------------------------------------------------------ */ diff --git a/src/MOLECULE/bond_harmonic.h b/src/MOLECULE/bond_harmonic.h index 443290aec1..75c1db2556 100644 --- a/src/MOLECULE/bond_harmonic.h +++ b/src/MOLECULE/bond_harmonic.h @@ -35,6 +35,7 @@ class BondHarmonic : public Bond { void read_restart(FILE *) override; void write_data(FILE *) override; double single(int, double, int, int, double &) override; + void born_matrix(int, double, int, int, double &, double &) override; void *extract(const char *, int &) override; protected: diff --git a/src/angle.cpp b/src/angle.cpp index a80129a087..6fdb9307cb 100644 --- a/src/angle.cpp +++ b/src/angle.cpp @@ -34,6 +34,7 @@ Angle::Angle(LAMMPS *_lmp) : Pointers(_lmp) allocated = 0; suffix_flag = Suffix::NONE; + born_matrix_enable = 0; maxeatom = maxvatom = maxcvatom = 0; eatom = nullptr; diff --git a/src/angle.h b/src/angle.h index d134f88aa2..ec10226b69 100644 --- a/src/angle.h +++ b/src/angle.h @@ -25,7 +25,8 @@ class Angle : protected Pointers { public: int allocated; int *setflag; - int writedata; // 1 if writes coeffs to data file + int writedata; // 1 if writes coeffs to data file + int born_matrix_enable; double energy; // accumulated energies double virial[6]; // accumulated virial: xx,yy,zz,xy,xz,yz double *eatom, **vatom; // accumulated per-atom energy/virial @@ -59,6 +60,12 @@ class Angle : protected Pointers { virtual void read_restart_settings(FILE *){}; virtual void write_data(FILE *) {} virtual double single(int, int, int, int) = 0; + virtual void born_matrix(int /*atype*/, int /*at1*/, int /*at2*/, int /*at3*/, double &du, + double &du2) + { + du = 0.0; + du2 = 0.0; + } virtual double memory_usage(); virtual void *extract(const char *, int &) { return nullptr; } void reinit(); diff --git a/src/atom.cpp b/src/atom.cpp index 69a3b21cf5..8283dbfee6 100644 --- a/src/atom.cpp +++ b/src/atom.cpp @@ -300,20 +300,20 @@ Atom::~Atom() // delete custom atom arrays for (int i = 0; i < nivector; i++) { - delete [] ivname[i]; + delete[] ivname[i]; memory->destroy(ivector[i]); } for (int i = 0; i < ndvector; i++) { - delete [] dvname[i]; + delete[] dvname[i]; if (dvector) // (needed for Kokkos) memory->destroy(dvector[i]); } for (int i = 0; i < niarray; i++) { - delete [] ianame[i]; + delete[] ianame[i]; memory->destroy(iarray[i]); } for (int i = 0; i < ndarray; i++) { - delete [] daname[i]; + delete[] daname[i]; memory->destroy(darray[i]); } @@ -335,8 +335,8 @@ Atom::~Atom() // delete per-type arrays - delete [] mass; - delete [] mass_setflag; + delete[] mass; + delete[] mass_setflag; // delete extra arrays @@ -730,7 +730,7 @@ void Atom::init() if (firstgroupname) { firstgroup = group->find(firstgroupname); if (firstgroup < 0) - error->all(FLERR,"Could not find atom_modify first group ID"); + error->all(FLERR,"Could not find atom_modify first group ID {}", firstgroupname); } else firstgroup = -1; // init AtomVec @@ -773,30 +773,30 @@ AtomVec *Atom::style_match(const char *style) void Atom::modify_params(int narg, char **arg) { - if (narg == 0) error->all(FLERR,"Illegal atom_modify command"); + if (narg == 0) utils::missing_cmd_args(FLERR, "atom_modify", error); int iarg = 0; while (iarg < narg) { if (strcmp(arg[iarg],"id") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal atom_modify command"); + if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "atom_modify id", error); if (domain->box_exist) error->all(FLERR,"Atom_modify id command after simulation box is defined"); tag_enable = utils::logical(FLERR,arg[iarg+1],false,lmp); iarg += 2; } else if (strcmp(arg[iarg],"map") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal atom_modify command"); + if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "atom_modify map", error); if (domain->box_exist) error->all(FLERR,"Atom_modify map command after simulation box is defined"); if (strcmp(arg[iarg+1],"array") == 0) map_user = 1; else if (strcmp(arg[iarg+1],"hash") == 0) map_user = 2; else if (strcmp(arg[iarg+1],"yes") == 0) map_user = 3; - else error->all(FLERR,"Illegal atom_modify command"); + else error->all(FLERR,"Illegal atom_modify map command argument {}", arg[iarg+1]); map_style = map_user; iarg += 2; } else if (strcmp(arg[iarg],"first") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal atom_modify command"); + if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "atom_modify first", error); if (strcmp(arg[iarg+1],"all") == 0) { - delete [] firstgroupname; + delete[] firstgroupname; firstgroupname = nullptr; } else { firstgroupname = utils::strdup(arg[iarg+1]); @@ -804,16 +804,15 @@ void Atom::modify_params(int narg, char **arg) } iarg += 2; } else if (strcmp(arg[iarg],"sort") == 0) { - if (iarg+3 > narg) error->all(FLERR,"Illegal atom_modify command"); + if (iarg+3 > narg) utils::missing_cmd_args(FLERR, "atom_modify sort", error); sortfreq = utils::inumeric(FLERR,arg[iarg+1],false,lmp); userbinsize = utils::numeric(FLERR,arg[iarg+2],false,lmp); - if (sortfreq < 0 || userbinsize < 0.0) - error->all(FLERR,"Illegal atom_modify command"); - if (sortfreq >= 0 && firstgroupname) - error->all(FLERR,"Atom_modify sort and first options " - "cannot be used together"); + if (sortfreq < 0) error->all(FLERR,"Illegal atom_modify sort frequency {}", sortfreq); + if (userbinsize < 0.0) error->all(FLERR,"Illegal atom_modify sort bin size {}", userbinsize); + if ((sortfreq >= 0) && firstgroupname) + error->all(FLERR,"Atom_modify sort and first options cannot be used together"); iarg += 3; - } else error->all(FLERR,"Illegal atom_modify command"); + } else error->all(FLERR,"Illegal atom_modify command argument: {}", arg[iarg]); } } @@ -884,7 +883,7 @@ void Atom::tag_extend() bigint notag_total; MPI_Allreduce(¬ag,¬ag_total,1,MPI_LMP_BIGINT,MPI_SUM,world); if (notag_total >= MAXTAGINT) - error->all(FLERR,"New atom IDs exceed maximum allowed ID"); + error->all(FLERR,"New atom IDs exceed maximum allowed ID {}", MAXTAGINT); bigint notag_sum; MPI_Scan(¬ag,¬ag_sum,1,MPI_LMP_BIGINT,MPI_SUM,world); @@ -1016,12 +1015,13 @@ void Atom::data_atoms(int n, char *buf, tagint id_offset, tagint mol_offset, // use the first line to detect and validate the number of words/tokens per line next = strchr(buf,'\n'); + if (!next) error->all(FLERR, "Missing data in Atoms section of data file"); *next = '\0'; int nwords = utils::trim_and_count_words(buf); *next = '\n'; - if (nwords != avec->size_data_atom && nwords != avec->size_data_atom + 3) - error->all(FLERR,"Incorrect atom format in data file"); + if ((nwords != avec->size_data_atom) && (nwords != avec->size_data_atom + 3)) + error->all(FLERR,"Incorrect atom format in data file: {}", utils::trim(buf)); // set bounds for my proc // if periodic and I am lo/hi proc, adjust bounds by EPSILON @@ -1093,54 +1093,57 @@ void Atom::data_atoms(int n, char *buf, tagint id_offset, tagint mol_offset, for (int i = 0; i < n; i++) { next = strchr(buf,'\n'); + if (!next) error->all(FLERR, "Missing data in Atoms section of data file"); *next = '\0'; auto values = Tokenizer(utils::trim_comment(buf)).as_vector(); - if ((int)values.size() != nwords) + if (values.size() == 0) { + // skip over empty or comment lines + } else if ((int)values.size() != nwords) { error->all(FLERR, "Incorrect atom format in data file: {}", utils::trim(buf)); - - int imx = 0, imy = 0, imz = 0; - if (imageflag) { - imx = utils::inumeric(FLERR,values[iptr],false,lmp); - imy = utils::inumeric(FLERR,values[iptr+1],false,lmp); - imz = utils::inumeric(FLERR,values[iptr+2],false,lmp); - if ((domain->dimension == 2) && (imz != 0)) - error->all(FLERR,"Z-direction image flag must be 0 for 2d-systems"); - if ((!domain->xperiodic) && (imx != 0)) { reset_image_flag[0] = true; imx = 0; } - if ((!domain->yperiodic) && (imy != 0)) { reset_image_flag[1] = true; imy = 0; } - if ((!domain->zperiodic) && (imz != 0)) { reset_image_flag[2] = true; imz = 0; } - } - imagedata = ((imageint) (imx + IMGMAX) & IMGMASK) | + } else { + int imx = 0, imy = 0, imz = 0; + if (imageflag) { + imx = utils::inumeric(FLERR,values[iptr],false,lmp); + imy = utils::inumeric(FLERR,values[iptr+1],false,lmp); + imz = utils::inumeric(FLERR,values[iptr+2],false,lmp); + if ((domain->dimension == 2) && (imz != 0)) + error->all(FLERR,"Z-direction image flag must be 0 for 2d-systems"); + if ((!domain->xperiodic) && (imx != 0)) { reset_image_flag[0] = true; imx = 0; } + if ((!domain->yperiodic) && (imy != 0)) { reset_image_flag[1] = true; imy = 0; } + if ((!domain->zperiodic) && (imz != 0)) { reset_image_flag[2] = true; imz = 0; } + } + imagedata = ((imageint) (imx + IMGMAX) & IMGMASK) | (((imageint) (imy + IMGMAX) & IMGMASK) << IMGBITS) | (((imageint) (imz + IMGMAX) & IMGMASK) << IMG2BITS); - xdata[0] = utils::numeric(FLERR,values[xptr],false,lmp); - xdata[1] = utils::numeric(FLERR,values[xptr+1],false,lmp); - xdata[2] = utils::numeric(FLERR,values[xptr+2],false,lmp); - if (shiftflag) { - xdata[0] += shift[0]; - xdata[1] += shift[1]; - xdata[2] += shift[2]; - } + xdata[0] = utils::numeric(FLERR,values[xptr],false,lmp); + xdata[1] = utils::numeric(FLERR,values[xptr+1],false,lmp); + xdata[2] = utils::numeric(FLERR,values[xptr+2],false,lmp); + if (shiftflag) { + xdata[0] += shift[0]; + xdata[1] += shift[1]; + xdata[2] += shift[2]; + } - domain->remap(xdata,imagedata); - if (triclinic) { - domain->x2lamda(xdata,lamda); - coord = lamda; - } else coord = xdata; + domain->remap(xdata,imagedata); + if (triclinic) { + domain->x2lamda(xdata,lamda); + coord = lamda; + } else coord = xdata; - if (coord[0] >= sublo[0] && coord[0] < subhi[0] && - coord[1] >= sublo[1] && coord[1] < subhi[1] && - coord[2] >= sublo[2] && coord[2] < subhi[2]) { - avec->data_atom(xdata,imagedata,values); - if (id_offset) tag[nlocal-1] += id_offset; - if (mol_offset) molecule[nlocal-1] += mol_offset; - if (type_offset) { - type[nlocal-1] += type_offset; - if (type[nlocal-1] > ntypes) - error->one(FLERR,"Invalid atom type in Atoms section of data file"); + if (coord[0] >= sublo[0] && coord[0] < subhi[0] && + coord[1] >= sublo[1] && coord[1] < subhi[1] && + coord[2] >= sublo[2] && coord[2] < subhi[2]) { + avec->data_atom(xdata,imagedata,values); + if (id_offset) tag[nlocal-1] += id_offset; + if (mol_offset) molecule[nlocal-1] += mol_offset; + if (type_offset) { + type[nlocal-1] += type_offset; + if (type[nlocal-1] > ntypes) + error->one(FLERR,"Invalid atom type in Atoms section of data file"); + } } } - buf = next + 1; } } @@ -1156,30 +1159,25 @@ void Atom::data_vels(int n, char *buf, tagint id_offset) int m; char *next; - next = strchr(buf,'\n'); - *next = '\0'; - int nwords = utils::trim_and_count_words(buf); - *next = '\n'; - - if (nwords != avec->size_data_vel) - error->all(FLERR,"Incorrect velocity format in data file"); - // loop over lines of atom velocities // tokenize the line into values // if I own atom tag, unpack its values for (int i = 0; i < n; i++) { next = strchr(buf,'\n'); + if (!next) error->all(FLERR, "Missing data in Velocities section of data file"); *next = '\0'; auto values = Tokenizer(utils::trim_comment(buf)).as_vector(); - if ((int)values.size() != nwords) - error->all(FLERR, "Incorrect atom format in data file: {}", utils::trim(buf)); - - tagint tagdata = utils::tnumeric(FLERR,values[0],false,lmp) + id_offset; - if (tagdata <= 0 || tagdata > map_tag_max) - error->one(FLERR,"Invalid atom ID in Velocities section of data file"); - if ((m = map(tagdata)) >= 0) avec->data_vel(m,values); - + if (values.size() == 0) { + // skip over empty or comment lines + } else if ((int)values.size() != avec->size_data_vel) { + error->all(FLERR, "Incorrect velocity format in data file: {}", utils::trim(buf)); + } else { + tagint tagdata = utils::tnumeric(FLERR,values[0],false,lmp) + id_offset; + if (tagdata <= 0 || tagdata > map_tag_max) + error->one(FLERR,"Invalid atom ID {} in Velocities section of data file: {}", tagdata, buf); + if ((m = map(tagdata)) >= 0) avec->data_vel(m,values); + } buf = next + 1; } } @@ -1202,47 +1200,51 @@ void Atom::data_bonds(int n, char *buf, int *count, tagint id_offset, for (int i = 0; i < n; i++) { next = strchr(buf,'\n'); + if (!next) error->all(FLERR, "Missing data in Bonds section of data file"); *next = '\0'; - try { - ValueTokenizer values(utils::trim_comment(buf)); - values.next_int(); - itype = values.next_int(); - atom1 = values.next_tagint(); - atom2 = values.next_tagint(); - if (values.has_next()) throw TokenizerException("Too many tokens",""); - } catch (TokenizerException &e) { - error->one(FLERR,"{} in {}: {}", e.what(), location, utils::trim(buf)); - } - if (id_offset) { - atom1 += id_offset; - atom2 += id_offset; - } - itype += type_offset; - - if ((atom1 <= 0) || (atom1 > map_tag_max) || - (atom2 <= 0) || (atom2 > map_tag_max) || (atom1 == atom2)) - error->one(FLERR,"Invalid atom ID in {}: {}", location, utils::trim(buf)); - if (itype <= 0 || itype > nbondtypes) - error->one(FLERR,"Invalid bond type in {}: {}", location, utils::trim(buf)); - if ((m = map(atom1)) >= 0) { - if (count) count[m]++; - else { - bond_type[m][num_bond[m]] = itype; - bond_atom[m][num_bond[m]] = atom2; - num_bond[m]++; - avec->data_bonds_post(m, num_bond[m], atom1, atom2, id_offset); + ValueTokenizer values(utils::trim_comment(buf)); + // skip over empty of comment lines + if (values.has_next()) { + try { + values.next_int(); + itype = values.next_int(); + atom1 = values.next_tagint(); + atom2 = values.next_tagint(); + if (values.has_next()) throw TokenizerException("Too many tokens",""); + } catch (TokenizerException &e) { + error->one(FLERR,"{} in {}: {}", e.what(), location, utils::trim(buf)); } - } - if (newton_bond == 0) { - if ((m = map(atom2)) >= 0) { + if (id_offset) { + atom1 += id_offset; + atom2 += id_offset; + } + itype += type_offset; + + if ((atom1 <= 0) || (atom1 > map_tag_max) || + (atom2 <= 0) || (atom2 > map_tag_max) || (atom1 == atom2)) + error->one(FLERR,"Invalid atom ID in {}: {}", location, utils::trim(buf)); + if (itype <= 0 || itype > nbondtypes) + error->one(FLERR,"Invalid bond type in {}: {}", location, utils::trim(buf)); + if ((m = map(atom1)) >= 0) { if (count) count[m]++; else { bond_type[m][num_bond[m]] = itype; - bond_atom[m][num_bond[m]] = atom1; + bond_atom[m][num_bond[m]] = atom2; num_bond[m]++; avec->data_bonds_post(m, num_bond[m], atom1, atom2, id_offset); } } + if (newton_bond == 0) { + if ((m = map(atom2)) >= 0) { + if (count) count[m]++; + else { + bond_type[m][num_bond[m]] = itype; + bond_atom[m][num_bond[m]] = atom1; + num_bond[m]++; + avec->data_bonds_post(m, num_bond[m], atom1, atom2, id_offset); + } + } + } } buf = next + 1; } @@ -1266,44 +1268,36 @@ void Atom::data_angles(int n, char *buf, int *count, tagint id_offset, for (int i = 0; i < n; i++) { next = strchr(buf,'\n'); + if (!next) error->all(FLERR, "Missing data in Angles section of data file"); *next = '\0'; - try { - ValueTokenizer values(utils::trim_comment(buf)); - values.next_int(); - itype = values.next_int(); - atom1 = values.next_tagint(); - atom2 = values.next_tagint(); - atom3 = values.next_tagint(); - if (values.has_next()) throw TokenizerException("Too many tokens",""); - } catch (TokenizerException &e) { - error->one(FLERR,"{} in {}: {}", e.what(), location, utils::trim(buf)); - } - if (id_offset) { - atom1 += id_offset; - atom2 += id_offset; - atom3 += id_offset; - } - itype += type_offset; - - if ((atom1 <= 0) || (atom1 > map_tag_max) || - (atom2 <= 0) || (atom2 > map_tag_max) || - (atom3 <= 0) || (atom3 > map_tag_max) || - (atom1 == atom2) || (atom1 == atom3) || (atom2 == atom3)) - error->one(FLERR,"Invalid atom ID in {}: {}", location, utils::trim(buf)); - if (itype <= 0 || itype > nangletypes) - error->one(FLERR,"Invalid angle type in {}: {}", location, utils::trim(buf)); - if ((m = map(atom2)) >= 0) { - if (count) count[m]++; - else { - angle_type[m][num_angle[m]] = itype; - angle_atom1[m][num_angle[m]] = atom1; - angle_atom2[m][num_angle[m]] = atom2; - angle_atom3[m][num_angle[m]] = atom3; - num_angle[m]++; + ValueTokenizer values(utils::trim_comment(buf)); + // skip over empty or comment lines + if (values.has_next()) { + try { + values.next_int(); + itype = values.next_int(); + atom1 = values.next_tagint(); + atom2 = values.next_tagint(); + atom3 = values.next_tagint(); + if (values.has_next()) throw TokenizerException("Too many tokens",""); + } catch (TokenizerException &e) { + error->one(FLERR,"{} in {}: {}", e.what(), location, utils::trim(buf)); } - } - if (newton_bond == 0) { - if ((m = map(atom1)) >= 0) { + if (id_offset) { + atom1 += id_offset; + atom2 += id_offset; + atom3 += id_offset; + } + itype += type_offset; + + if ((atom1 <= 0) || (atom1 > map_tag_max) || + (atom2 <= 0) || (atom2 > map_tag_max) || + (atom3 <= 0) || (atom3 > map_tag_max) || + (atom1 == atom2) || (atom1 == atom3) || (atom2 == atom3)) + error->one(FLERR,"Invalid atom ID in {}: {}", location, utils::trim(buf)); + if (itype <= 0 || itype > nangletypes) + error->one(FLERR,"Invalid angle type in {}: {}", location, utils::trim(buf)); + if ((m = map(atom2)) >= 0) { if (count) count[m]++; else { angle_type[m][num_angle[m]] = itype; @@ -1313,14 +1307,26 @@ void Atom::data_angles(int n, char *buf, int *count, tagint id_offset, num_angle[m]++; } } - if ((m = map(atom3)) >= 0) { - if (count) count[m]++; - else { - angle_type[m][num_angle[m]] = itype; - angle_atom1[m][num_angle[m]] = atom1; - angle_atom2[m][num_angle[m]] = atom2; - angle_atom3[m][num_angle[m]] = atom3; - num_angle[m]++; + if (newton_bond == 0) { + if ((m = map(atom1)) >= 0) { + if (count) count[m]++; + else { + angle_type[m][num_angle[m]] = itype; + angle_atom1[m][num_angle[m]] = atom1; + angle_atom2[m][num_angle[m]] = atom2; + angle_atom3[m][num_angle[m]] = atom3; + num_angle[m]++; + } + } + if ((m = map(atom3)) >= 0) { + if (count) count[m]++; + else { + angle_type[m][num_angle[m]] = itype; + angle_atom1[m][num_angle[m]] = atom1; + angle_atom2[m][num_angle[m]] = atom2; + angle_atom3[m][num_angle[m]] = atom3; + num_angle[m]++; + } } } } @@ -1346,49 +1352,40 @@ void Atom::data_dihedrals(int n, char *buf, int *count, tagint id_offset, for (int i = 0; i < n; i++) { next = strchr(buf,'\n'); + if (!next) error->all(FLERR, "Missing data in Dihedrals section of data file"); *next = '\0'; - try { - ValueTokenizer values(utils::trim_comment(buf)); - values.next_int(); - itype = values.next_int(); - atom1 = values.next_tagint(); - atom2 = values.next_tagint(); - atom3 = values.next_tagint(); - atom4 = values.next_tagint(); - if (values.has_next()) throw TokenizerException("Too many tokens",""); - } catch (TokenizerException &e) { - error->one(FLERR,"{} in {}: {}", e.what(), location, utils::trim(buf)); - } - if (id_offset) { - atom1 += id_offset; - atom2 += id_offset; - atom3 += id_offset; - atom4 += id_offset; - } - itype += type_offset; + ValueTokenizer values(utils::trim_comment(buf)); + // skip over empty or comment lines + if (values.has_next()) { + try { + values.next_int(); + itype = values.next_int(); + atom1 = values.next_tagint(); + atom2 = values.next_tagint(); + atom3 = values.next_tagint(); + atom4 = values.next_tagint(); + if (values.has_next()) throw TokenizerException("Too many tokens",""); + } catch (TokenizerException &e) { + error->one(FLERR,"{} in {}: {}", e.what(), location, utils::trim(buf)); + } + if (id_offset) { + atom1 += id_offset; + atom2 += id_offset; + atom3 += id_offset; + atom4 += id_offset; + } + itype += type_offset; - if ((atom1 <= 0) || (atom1 > map_tag_max) || - (atom2 <= 0) || (atom2 > map_tag_max) || - (atom3 <= 0) || (atom3 > map_tag_max) || - (atom4 <= 0) || (atom4 > map_tag_max) || - (atom1 == atom2) || (atom1 == atom3) || (atom1 == atom4) || - (atom2 == atom3) || (atom2 == atom4) || (atom3 == atom4)) - error->one(FLERR, "Invalid atom ID in {}: {}", location, utils::trim(buf)); - if (itype <= 0 || itype > ndihedraltypes) - error->one(FLERR, "Invalid dihedral type in {}: {}", location, utils::trim(buf)); - if ((m = map(atom2)) >= 0) { - if (count) count[m]++; - else { - dihedral_type[m][num_dihedral[m]] = itype; - dihedral_atom1[m][num_dihedral[m]] = atom1; - dihedral_atom2[m][num_dihedral[m]] = atom2; - dihedral_atom3[m][num_dihedral[m]] = atom3; - dihedral_atom4[m][num_dihedral[m]] = atom4; - num_dihedral[m]++; - } - } - if (newton_bond == 0) { - if ((m = map(atom1)) >= 0) { + if ((atom1 <= 0) || (atom1 > map_tag_max) || + (atom2 <= 0) || (atom2 > map_tag_max) || + (atom3 <= 0) || (atom3 > map_tag_max) || + (atom4 <= 0) || (atom4 > map_tag_max) || + (atom1 == atom2) || (atom1 == atom3) || (atom1 == atom4) || + (atom2 == atom3) || (atom2 == atom4) || (atom3 == atom4)) + error->one(FLERR, "Invalid atom ID in {}: {}", location, utils::trim(buf)); + if (itype <= 0 || itype > ndihedraltypes) + error->one(FLERR, "Invalid dihedral type in {}: {}", location, utils::trim(buf)); + if ((m = map(atom2)) >= 0) { if (count) count[m]++; else { dihedral_type[m][num_dihedral[m]] = itype; @@ -1399,26 +1396,39 @@ void Atom::data_dihedrals(int n, char *buf, int *count, tagint id_offset, num_dihedral[m]++; } } - if ((m = map(atom3)) >= 0) { - if (count) count[m]++; - else { - dihedral_type[m][num_dihedral[m]] = itype; - dihedral_atom1[m][num_dihedral[m]] = atom1; - dihedral_atom2[m][num_dihedral[m]] = atom2; - dihedral_atom3[m][num_dihedral[m]] = atom3; - dihedral_atom4[m][num_dihedral[m]] = atom4; - num_dihedral[m]++; + if (newton_bond == 0) { + if ((m = map(atom1)) >= 0) { + if (count) count[m]++; + else { + dihedral_type[m][num_dihedral[m]] = itype; + dihedral_atom1[m][num_dihedral[m]] = atom1; + dihedral_atom2[m][num_dihedral[m]] = atom2; + dihedral_atom3[m][num_dihedral[m]] = atom3; + dihedral_atom4[m][num_dihedral[m]] = atom4; + num_dihedral[m]++; + } } - } - if ((m = map(atom4)) >= 0) { - if (count) count[m]++; - else { - dihedral_type[m][num_dihedral[m]] = itype; - dihedral_atom1[m][num_dihedral[m]] = atom1; - dihedral_atom2[m][num_dihedral[m]] = atom2; - dihedral_atom3[m][num_dihedral[m]] = atom3; - dihedral_atom4[m][num_dihedral[m]] = atom4; - num_dihedral[m]++; + if ((m = map(atom3)) >= 0) { + if (count) count[m]++; + else { + dihedral_type[m][num_dihedral[m]] = itype; + dihedral_atom1[m][num_dihedral[m]] = atom1; + dihedral_atom2[m][num_dihedral[m]] = atom2; + dihedral_atom3[m][num_dihedral[m]] = atom3; + dihedral_atom4[m][num_dihedral[m]] = atom4; + num_dihedral[m]++; + } + } + if ((m = map(atom4)) >= 0) { + if (count) count[m]++; + else { + dihedral_type[m][num_dihedral[m]] = itype; + dihedral_atom1[m][num_dihedral[m]] = atom1; + dihedral_atom2[m][num_dihedral[m]] = atom2; + dihedral_atom3[m][num_dihedral[m]] = atom3; + dihedral_atom4[m][num_dihedral[m]] = atom4; + num_dihedral[m]++; + } } } } @@ -1444,49 +1454,40 @@ void Atom::data_impropers(int n, char *buf, int *count, tagint id_offset, for (int i = 0; i < n; i++) { next = strchr(buf,'\n'); + if (!next) error->all(FLERR, "Missing data in Impropers section of data file"); *next = '\0'; - try { - ValueTokenizer values(utils::trim_comment(buf)); - values.next_int(); - itype = values.next_int(); - atom1 = values.next_tagint(); - atom2 = values.next_tagint(); - atom3 = values.next_tagint(); - atom4 = values.next_tagint(); - if (values.has_next()) throw TokenizerException("Too many tokens",""); - } catch (TokenizerException &e) { - error->one(FLERR,"{} in {}: {}", e.what(), location, utils::trim(buf)); - } - if (id_offset) { - atom1 += id_offset; - atom2 += id_offset; - atom3 += id_offset; - atom4 += id_offset; - } - itype += type_offset; + ValueTokenizer values(utils::trim_comment(buf)); + // skip over empty or comment lines + if (values.has_next()) { + try { + values.next_int(); + itype = values.next_int(); + atom1 = values.next_tagint(); + atom2 = values.next_tagint(); + atom3 = values.next_tagint(); + atom4 = values.next_tagint(); + if (values.has_next()) throw TokenizerException("Too many tokens",""); + } catch (TokenizerException &e) { + error->one(FLERR,"{} in {}: {}", e.what(), location, utils::trim(buf)); + } + if (id_offset) { + atom1 += id_offset; + atom2 += id_offset; + atom3 += id_offset; + atom4 += id_offset; + } + itype += type_offset; - if ((atom1 <= 0) || (atom1 > map_tag_max) || - (atom2 <= 0) || (atom2 > map_tag_max) || - (atom3 <= 0) || (atom3 > map_tag_max) || - (atom4 <= 0) || (atom4 > map_tag_max) || - (atom1 == atom2) || (atom1 == atom3) || (atom1 == atom4) || - (atom2 == atom3) || (atom2 == atom4) || (atom3 == atom4)) - error->one(FLERR, "Invalid atom ID in {}: {}", location, utils::trim(buf)); - if (itype <= 0 || itype > nimpropertypes) - error->one(FLERR, "Invalid improper type in {}: {}", location, utils::trim(buf)); - if ((m = map(atom2)) >= 0) { - if (count) count[m]++; - else { - improper_type[m][num_improper[m]] = itype; - improper_atom1[m][num_improper[m]] = atom1; - improper_atom2[m][num_improper[m]] = atom2; - improper_atom3[m][num_improper[m]] = atom3; - improper_atom4[m][num_improper[m]] = atom4; - num_improper[m]++; - } - } - if (newton_bond == 0) { - if ((m = map(atom1)) >= 0) { + if ((atom1 <= 0) || (atom1 > map_tag_max) || + (atom2 <= 0) || (atom2 > map_tag_max) || + (atom3 <= 0) || (atom3 > map_tag_max) || + (atom4 <= 0) || (atom4 > map_tag_max) || + (atom1 == atom2) || (atom1 == atom3) || (atom1 == atom4) || + (atom2 == atom3) || (atom2 == atom4) || (atom3 == atom4)) + error->one(FLERR, "Invalid atom ID in {}: {}", location, utils::trim(buf)); + if (itype <= 0 || itype > nimpropertypes) + error->one(FLERR, "Invalid improper type in {}: {}", location, utils::trim(buf)); + if ((m = map(atom2)) >= 0) { if (count) count[m]++; else { improper_type[m][num_improper[m]] = itype; @@ -1497,26 +1498,39 @@ void Atom::data_impropers(int n, char *buf, int *count, tagint id_offset, num_improper[m]++; } } - if ((m = map(atom3)) >= 0) { - if (count) count[m]++; - else { - improper_type[m][num_improper[m]] = itype; - improper_atom1[m][num_improper[m]] = atom1; - improper_atom2[m][num_improper[m]] = atom2; - improper_atom3[m][num_improper[m]] = atom3; - improper_atom4[m][num_improper[m]] = atom4; - num_improper[m]++; + if (newton_bond == 0) { + if ((m = map(atom1)) >= 0) { + if (count) count[m]++; + else { + improper_type[m][num_improper[m]] = itype; + improper_atom1[m][num_improper[m]] = atom1; + improper_atom2[m][num_improper[m]] = atom2; + improper_atom3[m][num_improper[m]] = atom3; + improper_atom4[m][num_improper[m]] = atom4; + num_improper[m]++; + } } - } - if ((m = map(atom4)) >= 0) { - if (count) count[m]++; - else { - improper_type[m][num_improper[m]] = itype; - improper_atom1[m][num_improper[m]] = atom1; - improper_atom2[m][num_improper[m]] = atom2; - improper_atom3[m][num_improper[m]] = atom3; - improper_atom4[m][num_improper[m]] = atom4; - num_improper[m]++; + if ((m = map(atom3)) >= 0) { + if (count) count[m]++; + else { + improper_type[m][num_improper[m]] = itype; + improper_atom1[m][num_improper[m]] = atom1; + improper_atom2[m][num_improper[m]] = atom2; + improper_atom3[m][num_improper[m]] = atom3; + improper_atom4[m][num_improper[m]] = atom4; + num_improper[m]++; + } + } + if ((m = map(atom4)) >= 0) { + if (count) count[m]++; + else { + improper_type[m][num_improper[m]] = itype; + improper_atom1[m][num_improper[m]] = atom1; + improper_atom2[m][num_improper[m]] = atom2; + improper_atom3[m][num_improper[m]] = atom3; + improper_atom4[m][num_improper[m]] = atom4; + num_improper[m]++; + } } } } @@ -1535,34 +1549,29 @@ void Atom::data_bonus(int n, char *buf, AtomVec *avec_bonus, tagint id_offset) int m; char *next; - next = strchr(buf,'\n'); - *next = '\0'; - int nwords = utils::trim_and_count_words(buf); - *next = '\n'; - - if (nwords != avec_bonus->size_data_bonus) - error->all(FLERR,"Incorrect bonus data format in data file"); - // loop over lines of bonus atom data // tokenize the line into values // if I own atom tag, unpack its values for (int i = 0; i < n; i++) { next = strchr(buf,'\n'); + if (!next) error->all(FLERR, "Missing data in Bonus section of data file"); *next = '\0'; auto values = Tokenizer(utils::trim_comment(buf)).as_vector(); - if ((int)values.size() != nwords) - error->all(FLERR, "Incorrect atom format in data file: {}", utils::trim(buf)); + if (values.size() == 0) { + // skip over empty or comment lines + } else if ((int)values.size() != avec_bonus->size_data_bonus) { + error->all(FLERR, "Incorrect bonus data format in data file: {}", utils::trim(buf)); + } else { + tagint tagdata = utils::tnumeric(FLERR,values[0],false,lmp) + id_offset; + if (tagdata <= 0 || tagdata > map_tag_max) + error->one(FLERR,"Invalid atom ID in Bonus section of data file"); - tagint tagdata = utils::tnumeric(FLERR,values[0],false,lmp) + id_offset; - if (tagdata <= 0 || tagdata > map_tag_max) - error->one(FLERR,"Invalid atom ID in Bonus section of data file"); - - // ok to call child's data_atom_bonus() method thru parent avec_bonus, - // since data_bonus() was called with child ptr, and method is virtual - - if ((m = map(tagdata)) >= 0) avec_bonus->data_atom_bonus(m,values); + // ok to call child's data_atom_bonus() method thru parent avec_bonus, + // since data_bonus() was called with child ptr, and method is virtual + if ((m = map(tagdata)) >= 0) avec_bonus->data_atom_bonus(m,values); + } buf = next + 1; } } @@ -1587,46 +1596,49 @@ void Atom::data_bodies(int n, char *buf, AtomVec *avec_body, tagint id_offset) for (int i = 0; i < n; i++) { char *next = strchr(buf,'\n'); + if (!next) error->all(FLERR, "Missing data in Bodies section of data file"); *next = '\0'; auto values = Tokenizer(utils::trim_comment(buf)).as_vector(); - tagint tagdata = utils::tnumeric(FLERR,values[0],false,lmp) + id_offset; - int ninteger = utils::inumeric(FLERR,values[1],false,lmp); - int ndouble = utils::inumeric(FLERR,values[2],false,lmp); + if (values.size()) { + tagint tagdata = utils::tnumeric(FLERR,values[0],false,lmp) + id_offset; + int ninteger = utils::inumeric(FLERR,values[1],false,lmp); + int ndouble = utils::inumeric(FLERR,values[2],false,lmp); - if (unique_tags->find(tagdata) == unique_tags->end()) - unique_tags->insert(tagdata); - else - error->one(FLERR,"Duplicate atom ID in Bodies section of data file"); + if (unique_tags->find(tagdata) == unique_tags->end()) + unique_tags->insert(tagdata); + else + error->one(FLERR,"Duplicate atom ID {} in Bodies section of data file", tagdata); - buf = next + 1; - int m = map(tagdata); - if (m >= 0) { - ivalues.resize(ninteger); - dvalues.resize(ndouble); + buf = next + 1; + int m = map(tagdata); + if (m >= 0) { + ivalues.resize(ninteger); + dvalues.resize(ndouble); - for (int j = 0; j < ninteger; j++) { - buf += strspn(buf," \t\n\r\f"); - buf[strcspn(buf," \t\n\r\f")] = '\0'; - ivalues[j] = utils::inumeric(FLERR,buf,false,lmp); - buf += strlen(buf)+1; - } + for (int j = 0; j < ninteger; j++) { + buf += strspn(buf," \t\n\r\f"); + buf[strcspn(buf," \t\n\r\f")] = '\0'; + ivalues[j] = utils::inumeric(FLERR,buf,false,lmp); + buf += strlen(buf)+1; + } - for (int j = 0; j < ndouble; j++) { - buf += strspn(buf," \t\n\r\f"); - buf[strcspn(buf," \t\n\r\f")] = '\0'; - dvalues[j] = utils::numeric(FLERR,buf,false,lmp); - buf += strlen(buf)+1; - } + for (int j = 0; j < ndouble; j++) { + buf += strspn(buf," \t\n\r\f"); + buf[strcspn(buf," \t\n\r\f")] = '\0'; + dvalues[j] = utils::numeric(FLERR,buf,false,lmp); + buf += strlen(buf)+1; + } - avec_body->data_body(m,ninteger,ndouble,ivalues.data(),dvalues.data()); + avec_body->data_body(m,ninteger,ndouble,ivalues.data(),dvalues.data()); - } else { - int nvalues = ninteger + ndouble; // number of values to skip - for (int j = 0; j < nvalues; j++) { - buf += strspn(buf," \t\n\r\f"); - buf[strcspn(buf," \t\n\r\f")] = '\0'; - buf += strlen(buf)+1; + } else { + int nvalues = ninteger + ndouble; // number of values to skip + for (int j = 0; j < nvalues; j++) { + buf += strspn(buf," \t\n\r\f"); + buf[strcspn(buf," \t\n\r\f")] = '\0'; + buf += strlen(buf)+1; + } } } buf += strspn(buf," \t\n\r\f"); @@ -1681,24 +1693,26 @@ void Atom::allocate_type_arrays() void Atom::set_mass(const char *file, int line, const char *str, int type_offset) { - if (mass == nullptr) error->all(file,line,"Cannot set mass for this atom style"); + if (mass == nullptr) error->all(file,line,"Cannot set mass for atom style {}", atom_style); int itype; double mass_one; - try { - ValueTokenizer values(utils::trim_comment(str)); - itype = values.next_int() + type_offset; - mass_one = values.next_double(); - if (values.has_next()) throw TokenizerException("Too many tokens", ""); + ValueTokenizer values(utils::trim_comment(str)); + if (values.has_next()) { + try { + itype = values.next_int() + type_offset; + mass_one = values.next_double(); + if (values.has_next()) throw TokenizerException("Too many tokens", ""); - if (itype < 1 || itype > ntypes) throw TokenizerException("Invalid atom type", ""); - if (mass_one <= 0.0) throw TokenizerException("Invalid mass value", ""); - } catch (TokenizerException &e) { - error->all(file,line,"{} in Masses section of data file: {}", e.what(), utils::trim(str)); + if (itype < 1 || itype > ntypes) throw TokenizerException("Invalid atom type", ""); + if (mass_one <= 0.0) throw TokenizerException("Invalid mass value", ""); + } catch (TokenizerException &e) { + error->all(file,line,"{} in Masses section of data file: {}", e.what(), utils::trim(str)); + } + + mass[itype] = mass_one; + mass_setflag[itype] = 1; } - - mass[itype] = mass_one; - mass_setflag[itype] = 1; } /* ---------------------------------------------------------------------- @@ -1708,14 +1722,13 @@ void Atom::set_mass(const char *file, int line, const char *str, int type_offset void Atom::set_mass(const char *file, int line, int itype, double value) { - if (mass == nullptr) error->all(file,line,"Cannot set mass for this atom style"); + if (mass == nullptr) error->all(file,line,"Cannot set mass for atom style {}", atom_style); if (itype < 1 || itype > ntypes) - error->all(file,line,"Invalid type for mass set"); + error->all(file,line,"Invalid type {} for atom mass {}", itype, value); + if (value <= 0.0) error->all(file,line,"Invalid atom mass value {}", value); mass[itype] = value; mass_setflag[itype] = 1; - - if (mass[itype] <= 0.0) error->all(file,line,"Invalid mass value"); } /* ---------------------------------------------------------------------- @@ -1725,17 +1738,19 @@ void Atom::set_mass(const char *file, int line, int itype, double value) void Atom::set_mass(const char *file, int line, int /*narg*/, char **arg) { - if (mass == nullptr) error->all(file,line,"Cannot set mass for this atom style"); + if (mass == nullptr) error->all(file,line,"Cannot set atom mass for atom style", atom_style); int lo,hi; utils::bounds(file,line,arg[0],1,ntypes,lo,hi,error); - if (lo < 1 || hi > ntypes) error->all(file,line,"Invalid type for mass set"); + if ((lo < 1) || (hi > ntypes)) + error->all(file,line,"Invalid type {} for atom mass {}", arg[1]); + + const double value = utils::numeric(FLERR,arg[1],false,lmp); + if (value <= 0.0) error->all(file,line,"Invalid atom mass value {}", value); for (int itype = lo; itype <= hi; itype++) { - mass[itype] = utils::numeric(FLERR,arg[1],false,lmp); + mass[itype] = value; mass_setflag[itype] = 1; - - if (mass[itype] <= 0.0) error->all(file,line,"Invalid mass value"); } } @@ -1762,7 +1777,7 @@ void Atom::check_mass(const char *file, int line) if (rmass_flag) return; for (int itype = 1; itype <= ntypes; itype++) if (mass_setflag[itype] == 0) - error->all(file,line,"Not all per-type masses are set"); + error->all(file,line,"Not all per-type masses are set. Type {} is missing.", itype); } /* ---------------------------------------------------------------------- @@ -1795,8 +1810,7 @@ int Atom::radius_consistency(int itype, double &rad) also return the 3 shape params for itype ------------------------------------------------------------------------- */ -int Atom::shape_consistency(int itype, - double &shapex, double &shapey, double &shapez) +int Atom::shape_consistency(int itype, double &shapex, double &shapey, double &shapez) { double zero[3] = {0.0, 0.0, 0.0}; double one[3] = {-1.0, -1.0, -1.0}; @@ -1815,7 +1829,7 @@ int Atom::shape_consistency(int itype, one[0] = shape[0]; one[1] = shape[1]; one[2] = shape[2]; - } else if (one[0] != shape[0] || one[1] != shape[1] || one[2] != shape[2]) + } else if ((one[0] != shape[0]) || (one[1] != shape[1]) || (one[2] != shape[2])) flag = 1; } @@ -1837,10 +1851,10 @@ int Atom::shape_consistency(int itype, void Atom::add_molecule(int narg, char **arg) { - if (narg < 1) error->all(FLERR,"Illegal molecule command"); + if (narg < 1) utils::missing_cmd_args(FLERR, "molecule", error); if (find_molecule(arg[0]) >= 0) - error->all(FLERR,"Reuse of molecule template ID"); + error->all(FLERR,"Reuse of molecule template ID {}", arg[0]); // 1st molecule in set stores nset = # of mols, others store nset = 0 // ifile = count of molecules in set @@ -2110,8 +2124,7 @@ void Atom::setup_sort_bins() if ((binsize == 0.0) && (sortfreq > 0)) { sortfreq = 0; if (comm->me == 0) - error->warning(FLERR,"No pairwise cutoff or binsize set. " - "Atom sorting therefore disabled."); + error->warning(FLERR,"No pairwise cutoff or binsize set. Atom sorting therefore disabled."); return; } @@ -2213,8 +2226,7 @@ void Atom::setup_sort_bins() } #endif - if (1.0*nbinx*nbiny*nbinz > INT_MAX) - error->one(FLERR,"Too many atom sorting bins"); + if (1.0*nbinx*nbiny*nbinz > INT_MAX) error->one(FLERR,"Too many atom sorting bins"); nbins = nbinx*nbiny*nbinz; @@ -2459,25 +2471,25 @@ void Atom::remove_custom(int index, int flag, int cols) if (flag == 0 && cols == 0) { memory->destroy(ivector[index]); ivector[index] = nullptr; - delete [] ivname[index]; + delete[] ivname[index]; ivname[index] = nullptr; } else if (flag == 1 && cols == 0) { memory->destroy(dvector[index]); dvector[index] = nullptr; - delete [] dvname[index]; + delete[] dvname[index]; dvname[index] = nullptr; } else if (flag == 0 && cols) { memory->destroy(iarray[index]); iarray[index] = nullptr; - delete [] ianame[index]; + delete[] ianame[index]; ianame[index] = nullptr; } else if (flag == 1 && cols) { memory->destroy(darray[index]); darray[index] = nullptr; - delete [] daname[index]; + delete[] daname[index]; daname[index] = nullptr; } } diff --git a/src/bond.cpp b/src/bond.cpp index 2d747ca503..18cdaf8c6d 100644 --- a/src/bond.cpp +++ b/src/bond.cpp @@ -49,6 +49,7 @@ Bond::Bond(LAMMPS *_lmp) : Pointers(_lmp) allocated = 0; suffix_flag = Suffix::NONE; + born_matrix_enable = 0; partial_flag = 0; maxeatom = maxvatom = 0; diff --git a/src/bond.h b/src/bond.h index 75110f228c..3c5cb0f2d2 100644 --- a/src/bond.h +++ b/src/bond.h @@ -33,6 +33,8 @@ class Bond : protected Pointers { double virial[6]; // accumulated virial: xx,yy,zz,xy,xz,yz double *eatom, **vatom; // accumulated per-atom energy/virial + int born_matrix_enable; + int comm_forward; // size of forward communication (0 if none) int comm_reverse; // size of reverse communication (0 if none) int comm_reverse_off; // size of reverse comm even if newton off @@ -69,6 +71,13 @@ class Bond : protected Pointers { virtual int pack_reverse_comm(int, int, double *) { return 0; } virtual void unpack_reverse_comm(int, int *, double *) {} + virtual void born_matrix(int /*btype*/, double /*rsq*/, int /*at1*/, int /*at2*/, double &du, + double &du2) + { + du = 0.0; + du2 = 0.0; + } + void write_file(int, char **); protected: diff --git a/src/compute.h b/src/compute.h index 3d6323b3d1..6d0e4bd0f1 100644 --- a/src/compute.h +++ b/src/compute.h @@ -108,7 +108,7 @@ class Compute : protected Pointers { Compute(class LAMMPS *, int, char **); ~Compute() override; void modify_params(int, char **); - void reset_extra_dof(); + virtual void reset_extra_dof(); virtual void init() = 0; virtual void init_list(int, class NeighList *) {} diff --git a/src/compute_temp_profile.cpp b/src/compute_temp_profile.cpp index ff6bc36a08..4b8e77029b 100644 --- a/src/compute_temp_profile.cpp +++ b/src/compute_temp_profile.cpp @@ -119,6 +119,9 @@ ComputeTempProfile::ComputeTempProfile(LAMMPS *lmp, int narg, char **arg) : nbins = nbinx*nbiny*nbinz; if (nbins <= 0) error->all(FLERR,"Illegal compute temp/profile command"); + nstreaming = (xflag==0 ? 0 : 1) + (yflag==0 ? 0 : 1) + (zflag==0 ? 0 : 1); + reset_extra_dof(); + memory->create(vbin,nbins,ncount,"temp/profile:vbin"); memory->create(binave,nbins,ncount,"temp/profile:binave"); @@ -197,9 +200,10 @@ void ComputeTempProfile::dof_compute() natoms_temp = group->count(igroup); dof = domain->dimension * natoms_temp; - // subtract additional d*Nbins DOF, as in Evans and Morriss paper + // subtract additional Nbins DOF for each adjusted direction, + // as in Evans and Morriss paper - dof -= extra_dof + fix_dof + domain->dimension*nbins; + dof -= extra_dof + fix_dof + nstreaming*nbins; if (dof > 0) tfactor = force->mvv2e / (dof * force->boltz); else tfactor = 0.0; } @@ -334,14 +338,19 @@ void ComputeTempProfile::compute_array() MPI_Allreduce(tbin,tbinall,nbins,MPI_DOUBLE,MPI_SUM,world); - int nper = domain->dimension; + double totcount = 0.0; for (i = 0; i < nbins; i++) { array[i][0] = binave[i][ncount-1]; + totcount += array[i][0]; + } + double nper = domain->dimension - (extra_dof + fix_dof)/totcount; + double dofbin, tfactorbin; + for (i = 0; i < nbins; i++) { if (array[i][0] > 0.0) { - dof = nper*array[i][0] - extra_dof; - if (dof > 0) tfactor = force->mvv2e / (dof * force->boltz); - else tfactor = 0.0; - array[i][1] = tfactor*tbinall[i]; + dofbin = nper*array[i][0] - nstreaming; + if (dofbin > 0) tfactorbin = force->mvv2e / (dofbin * force->boltz); + else tfactorbin = 0.0; + array[i][1] = tfactorbin*tbinall[i]; } else array[i][1] = 0.0; } } @@ -576,6 +585,12 @@ void ComputeTempProfile::bin_assign() /* ---------------------------------------------------------------------- */ +void ComputeTempProfile::reset_extra_dof() { + extra_dof = domain->dimension - nstreaming; +} + +/* ---------------------------------------------------------------------- */ + double ComputeTempProfile::memory_usage() { double bytes = (double)maxatom * sizeof(int); diff --git a/src/compute_temp_profile.h b/src/compute_temp_profile.h index f9f2e42156..2d73781e3f 100644 --- a/src/compute_temp_profile.h +++ b/src/compute_temp_profile.h @@ -34,6 +34,7 @@ class ComputeTempProfile : public Compute { void compute_vector() override; void compute_array() override; + void reset_extra_dof() override; void remove_bias(int, double *) override; void remove_bias_thr(int, double *, double *) override; void remove_bias_all() override; @@ -47,6 +48,7 @@ class ComputeTempProfile : public Compute { int nbinx, nbiny, nbinz, nbins; int ivx, ivy, ivz; double tfactor; + double nstreaming; int box_change, triclinic; int *periodicity; diff --git a/src/dihedral.cpp b/src/dihedral.cpp index 5b4f923273..b82e1ed84e 100644 --- a/src/dihedral.cpp +++ b/src/dihedral.cpp @@ -35,6 +35,7 @@ Dihedral::Dihedral(LAMMPS *_lmp) : Pointers(_lmp) allocated = 0; suffix_flag = Suffix::NONE; + born_matrix_enable = 0; maxeatom = maxvatom = maxcvatom = 0; eatom = nullptr; diff --git a/src/dihedral.h b/src/dihedral.h index 7bb7eb2650..5b7e62e800 100644 --- a/src/dihedral.h +++ b/src/dihedral.h @@ -25,7 +25,8 @@ class Dihedral : protected Pointers { public: int allocated; int *setflag; - int writedata; // 1 if writes coeffs to data file + int writedata; // 1 if writes coeffs to data file + int born_matrix_enable; double energy; // accumulated energy double virial[6]; // accumulated virial: xx,yy,zz,xy,xz,yz double *eatom, **vatom; // accumulated per-atom energy/virial @@ -55,6 +56,12 @@ class Dihedral : protected Pointers { virtual void read_restart_settings(FILE *){}; virtual void write_data(FILE *) {} virtual double memory_usage(); + virtual void born_matrix(int /*dtype*/, int /*at1*/, int /*at2*/, int /*at3*/, int /*at4*/, + double &du, double &du2) + { + du = 0.0; + du2 = 0.0; + } protected: int suffix_flag; // suffix compatibility flag diff --git a/src/fix_bond_history.cpp b/src/fix_bond_history.cpp index 67489e6402..a0af563b3b 100644 --- a/src/fix_bond_history.cpp +++ b/src/fix_bond_history.cpp @@ -251,7 +251,7 @@ void FixBondHistory::post_neighbor() double FixBondHistory::memory_usage() { - return maxbond * ndata * sizeof(double); + return (double) maxbond * ndata * sizeof(double); } /* ---------------------------------------------------------------------- */ diff --git a/src/improper.cpp b/src/improper.cpp index 7ba39368b7..bb3c4d12b5 100644 --- a/src/improper.cpp +++ b/src/improper.cpp @@ -33,6 +33,7 @@ Improper::Improper(LAMMPS *_lmp) : Pointers(_lmp) allocated = 0; suffix_flag = Suffix::NONE; + born_matrix_enable = 0; maxeatom = maxvatom = maxcvatom = 0; eatom = nullptr; diff --git a/src/improper.h b/src/improper.h index dbe8bee3f1..dacdfc9797 100644 --- a/src/improper.h +++ b/src/improper.h @@ -25,7 +25,8 @@ class Improper : protected Pointers { public: int allocated; int *setflag; - int writedata; // 1 if writes coeffs to data file + int writedata; // 1 if writes coeffs to data file + int born_matrix_enable; double energy; // accumulated energies double virial[6]; // accumulated virial: xx,yy,zz,xy,xz,yz double *eatom, **vatom; // accumulated per-atom energy/virial @@ -55,6 +56,12 @@ class Improper : protected Pointers { virtual void read_restart_settings(FILE *){}; virtual void write_data(FILE *) {} virtual double memory_usage(); + virtual void born_matrix(int /*dtype*/, int /*at1*/, int /*at2*/, int /*at3*/, int /*at4*/, + double &du, double &du2) + { + du = 0.0; + du2 = 0.0; + } protected: int suffix_flag; // suffix compatibility flag diff --git a/src/pair.cpp b/src/pair.cpp index aee89db49e..a8168bd4d7 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 f4aa82a586..3930847ff5 100644 --- a/src/pair.h +++ b/src/pair.h @@ -52,6 +52,7 @@ class Pair : protected Pointers { int comm_reverse_off; // size of reverse comm even if newton off int single_enable; // 1 if single() routine exists + int born_matrix_enable; // 1 if born_matrix() routine exists int single_hessian_enable; // 1 if single_hessian() routine exists int restartinfo; // 1 if pair style writes restart info int respa_enable; // 1 if inner/middle/outer rRESPA routines @@ -169,6 +170,12 @@ 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 = du2 = 0.0; + } + virtual void settings(int, char **) = 0; virtual void coeff(int, char **) = 0; diff --git a/src/pair_hybrid.cpp b/src/pair_hybrid.cpp index 82dce4d5f2..d0a3a2b053 100644 --- a/src/pair_hybrid.cpp +++ b/src/pair_hybrid.cpp @@ -391,6 +391,7 @@ void PairHybrid::flags() // single_enable = 1 if all sub-styles are set // respa_enable = 1 if all sub-styles are set // manybody_flag = 1 if any sub-style is set + // born_matrix_enable = 1 if all sub-styles are set // no_virial_fdotr_compute = 1 if any sub-style is set // ghostneigh = 1 if any sub-style is set // ewaldflag, pppmflag, msmflag, dipoleflag, dispersionflag, tip4pflag = 1 @@ -401,11 +402,13 @@ void PairHybrid::flags() compute_flag = 0; respa_enable = 0; restartinfo = 0; + born_matrix_enable = 0; for (m = 0; m < nstyles; m++) { if (styles[m]->single_enable) ++single_enable; if (styles[m]->respa_enable) ++respa_enable; if (styles[m]->restartinfo) ++restartinfo; + if (styles[m]->born_matrix_enable) ++born_matrix_enable; if (styles[m]->manybody_flag) manybody_flag = 1; if (styles[m]->no_virial_fdotr_compute) no_virial_fdotr_compute = 1; if (styles[m]->ghostneigh) ghostneigh = 1; @@ -422,6 +425,7 @@ void PairHybrid::flags() single_enable = (single_enable == nstyles) ? 1 : 0; respa_enable = (respa_enable == nstyles) ? 1 : 0; restartinfo = (restartinfo == nstyles) ? 1 : 0; + born_matrix_enable = (born_matrix_enable == nstyles) ? 1 : 0; init_svector(); // set centroidstressflag for pair hybrid @@ -867,6 +871,40 @@ double PairHybrid::single(int i, int j, int itype, int jtype, return esum; } +/* ---------------------------------------------------------------------- + call sub-style to compute born matrix interaction + error if sub-style does not support born_matrix call + since overlay could have multiple sub-styles, sum results explicitly +------------------------------------------------------------------------- */ + +void PairHybrid::born_matrix(int i, int j, int itype, int jtype, double rsq, + double factor_coul, double factor_lj, + double &dupair, double &du2pair) +{ + if (nmap[itype][jtype] == 0) + error->one(FLERR,"Invoked pair born_matrix on pair style none"); + + double du, du2; + dupair = du2pair = 0.0; + + for (int m = 0; m < nmap[itype][jtype]; m++) { + if (rsq < styles[map[itype][jtype][m]]->cutsq[itype][jtype]) { + if (styles[map[itype][jtype][m]]->born_matrix_enable == 0) + error->one(FLERR,"Pair hybrid sub-style does not support born_matrix call"); + + if ((special_lj[map[itype][jtype][m]] != nullptr) || + (special_coul[map[itype][jtype][m]] != nullptr)) + error->one(FLERR,"Pair hybrid born_matrix calls do not support" + " per sub-style special bond values"); + + du = du2 = 0.0; + styles[map[itype][jtype][m]]->born_matrix(i,j,itype,jtype,rsq,factor_coul,factor_lj,du,du2); + dupair += du; + du2pair += du2; + } + } +} + /* ---------------------------------------------------------------------- copy Pair::svector data ------------------------------------------------------------------------- */ diff --git a/src/pair_hybrid.h b/src/pair_hybrid.h index ee56783700..675bd16043 100644 --- a/src/pair_hybrid.h +++ b/src/pair_hybrid.h @@ -49,6 +49,8 @@ class PairHybrid : public Pair { void write_restart(FILE *) override; void read_restart(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 modify_params(int narg, char **arg) override; double memory_usage() override; diff --git a/src/pair_hybrid_scaled.cpp b/src/pair_hybrid_scaled.cpp index 49ebaabe24..167fd09cf7 100644 --- a/src/pair_hybrid_scaled.cpp +++ b/src/pair_hybrid_scaled.cpp @@ -413,7 +413,7 @@ double PairHybridScaled::single(int i, int j, int itype, int jtype, double rsq, (special_coul[map[itype][jtype][m]] != nullptr)) error->one(FLERR, "Pair hybrid single() does not support per sub-style special_bond"); - scale = scaleval[map[itype][jtype][m]]; + double scale = scaleval[map[itype][jtype][m]]; esum += scale * pstyle->single(i, j, itype, jtype, rsq, factor_coul, factor_lj, fone); fforce += scale * fone; } @@ -423,6 +423,57 @@ double PairHybridScaled::single(int i, int j, int itype, int jtype, double rsq, return esum; } +/* ---------------------------------------------------------------------- + call sub-style to compute born matrix interaction + error if sub-style does not support born_matrix call + since overlay could have multiple sub-styles, sum results explicitly +------------------------------------------------------------------------- */ + +void PairHybridScaled::born_matrix(int i, int j, int itype, int jtype, double rsq, + double factor_coul, double factor_lj, double &dupair, + double &du2pair) +{ + if (nmap[itype][jtype] == 0) error->one(FLERR, "Invoked pair born_matrix on pair style none"); + + // update scale values from variables where needed + + const int nvars = scalevars.size(); + if (nvars > 0) { + double *vals = new double[nvars]; + for (int k = 0; k < nvars; ++k) { + int m = input->variable->find(scalevars[k].c_str()); + if (m < 0) + error->all(FLERR, "Variable '{}' not found when updating scale factors", scalevars[k]); + vals[k] = input->variable->compute_equal(m); + } + for (int k = 0; k < nstyles; ++k) { + if (scaleidx[k] >= 0) scaleval[k] = vals[scaleidx[k]]; + } + delete[] vals; + } + + double du, du2, scale; + dupair = du2pair = scale = 0.0; + + for (int m = 0; m < nmap[itype][jtype]; m++) { + auto pstyle = styles[map[itype][jtype][m]]; + if (rsq < pstyle->cutsq[itype][jtype]) { + if (pstyle->born_matrix_enable == 0) + error->one(FLERR, "Pair hybrid sub-style does not support born_matrix call"); + + if ((special_lj[map[itype][jtype][m]] != nullptr) || + (special_coul[map[itype][jtype][m]] != nullptr)) + error->one(FLERR, "Pair hybrid born_matrix() does not support per sub-style special_bond"); + + du = du2 = 0.0; + scale = scaleval[map[itype][jtype][m]]; + pstyle->born_matrix(i, j, itype, jtype, rsq, factor_coul, factor_lj, du, du2); + dupair += scale * du; + du2pair += scale * du2; + } + } +} + /* ---------------------------------------------------------------------- set coeffs for one or more type pairs ------------------------------------------------------------------------- */ diff --git a/src/pair_hybrid_scaled.h b/src/pair_hybrid_scaled.h index 9bb8901846..a7789c7164 100644 --- a/src/pair_hybrid_scaled.h +++ b/src/pair_hybrid_scaled.h @@ -38,6 +38,7 @@ class PairHybridScaled : public PairHybrid { void write_restart(FILE *) override; void read_restart(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 init_svector() override; void copy_svector(int, int) override; diff --git a/src/pair_lj_cut.cpp b/src/pair_lj_cut.cpp index 503d7a3b33..73d46b26ce 100644 --- a/src/pair_lj_cut.cpp +++ b/src/pair_lj_cut.cpp @@ -39,6 +39,7 @@ using namespace MathConst; PairLJCut::PairLJCut(LAMMPS *lmp) : Pair(lmp) { respa_enable = 1; + born_matrix_enable = 1; writedata = 1; } @@ -672,6 +673,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; diff --git a/src/utils.cpp b/src/utils.cpp index 0d0bc91227..3bc8c761c6 100644 --- a/src/utils.cpp +++ b/src/utils.cpp @@ -120,6 +120,12 @@ std::string utils::strfind(const std::string &text, const std::string &pattern) return ""; } +void utils::missing_cmd_args(const std::string &file, int line, const std::string &cmd, + Error *error) +{ + if (error) error->all(file, line, "Illegal {} command: missing argument(s)", cmd); +} + /* specialization for the case of just a single string argument */ void utils::logmesg(LAMMPS *lmp, const std::string &mesg) @@ -571,14 +577,14 @@ void utils::bounds(const char *file, int line, const std::string &str, error->all(file, line, fmt::format("Invalid range string: {}", str)); if (nlo < nmin) - error->all(file, line, fmt::format("Numeric index {} is out of bounds " - "({}-{})", nlo, nmin, nmax)); + error->all(file, line, fmt::format("Numeric index {} is out of bounds ({}-{})", + nlo, nmin, nmax)); else if (nhi > nmax) - error->all(file, line, fmt::format("Numeric index {} is out of bounds " - "({}-{})", nhi, nmin, nmax)); + error->all(file, line, fmt::format("Numeric index {} is out of bounds ({}-{})", + nhi, nmin, nmax)); else if (nlo > nhi) - error->all(file, line, fmt::format("Numeric index {} is out of bounds " - "({}-{})", nlo, nmin, nhi)); + error->all(file, line, fmt::format("Numeric index {} is out of bounds ({}-{})", + nlo, nmin, nhi)); } } diff --git a/src/utils.h b/src/utils.h index 4dbb2fb8bb..08b9faea3b 100644 --- a/src/utils.h +++ b/src/utils.h @@ -47,6 +47,18 @@ namespace utils { std::string strfind(const std::string &text, const std::string &pattern); + /*! Print error message about missing arguments for command + * + * This function simplifies the repetitive reporting missing arguments to a command. + * + * \param file name of source file for error message + * \param line line number in source file for error message + * \param cmd name of the failing command + * \param error pointer to Error class instance (for abort) or nullptr */ + + [[noreturn]] void missing_cmd_args(const std::string &file, int line, const std::string &cmd, + Error *error); + /* Internal function handling the argument list for logmesg(). */ void fmtargs_logmesg(LAMMPS *lmp, fmt::string_view format, fmt::format_args args);