Add HMA compute
This commit is contained in:
Andrew Schultz
@ -173,65 +173,43 @@ There are also additional accelerated compute styles included in the
|
||||
LAMMPS distribution for faster performance on CPUs, GPUs, and KNLs.
|
||||
The individual style names on the "Commands
|
||||
compute"_Commands_compute.html doc page are followed by one or more of
|
||||
(g,i,k,o,t) to indicate which accelerated styles exist.
|
||||
(g,i,k,o,t) to indicate which accerlerated styles exist.
|
||||
|
||||
"ackland/atom"_compute_ackland_atom.html -
|
||||
"aggregate/atom"_compute_cluster_atom.html - aggregate ID for each atom
|
||||
"angle"_compute_angle.html -
|
||||
"angle/local"_compute_angle_local.html -
|
||||
"angle/local"_compute_bond_local.html - theta and energy of each angle
|
||||
"angmom/chunk"_compute_angmom_chunk.html - angular momentum for each chunk
|
||||
"basal/atom"_compute_basal_atom.html -
|
||||
"body/local"_compute_body_local.html - attributes of body sub-particles
|
||||
"bond"_compute_bond.html - values computed by a bond style
|
||||
"bond/local"_compute_bond_local.html - distance and energy of each bond
|
||||
"centro/atom"_compute_centro_atom.html - centro-symmetry parameter for each atom
|
||||
"chunk/atom"_compute_chunk_atom.html - assign chunk IDs to each atom
|
||||
"chunk/spread/atom"_compute_chunk_spread_atom.html - spreads chunk values to each atom in chunk
|
||||
"cluster/atom"_compute_cluster_atom.html - cluster ID for each atom
|
||||
"cna/atom"_compute_cna_atom.html - common neighbor analysis (CNA) for each atom
|
||||
"cnp/atom"_compute_cnp_atom.html -
|
||||
"com"_compute_com.html - center-of-mass of group of atoms
|
||||
"com/chunk"_compute_com_chunk.html - center-of-mass for each chunk
|
||||
"contact/atom"_compute_contact_atom.html - contact count for each spherical particle
|
||||
"coord/atom"_compute_coord_atom.html - coordination number for each atom
|
||||
"damage/atom"_compute_damage_atom.html - Peridynamic damage for each atom
|
||||
"dihedral"_compute_dihedral.html -
|
||||
"dihedral/local"_compute_dihedral_local.html - angle of each dihedral
|
||||
"dilatation/atom"_compute_dilatation_atom.html - Peridynamic dilatation for each atom
|
||||
"dipole/chunk"_compute_dipole_chunk.html -
|
||||
"displace/atom"_compute_displace_atom.html - displacement of each atom
|
||||
"dpd"_compute_dpd.html -
|
||||
"dpd/atom"_compute_dpd_atom.html -
|
||||
"edpd/temp/atom"_compute_edpd_temp_atom.html -
|
||||
"entropy/atom"_compute_entropy_atom.html -
|
||||
"erotate/asphere"_compute_erotate_asphere.html - rotational energy of aspherical particles
|
||||
"erotate/rigid"_compute_erotate_rigid.html - rotational energy of rigid bodies
|
||||
"erotate/sphere"_compute_erotate_sphere.html - rotational energy of spherical particles
|
||||
"erotate/sphere/atom"_compute_erotate_sphere.html - rotational energy for each spherical particle
|
||||
"erotate/sphere/atom"_compute_erotate_sphere_atom.html -
|
||||
"event/displace"_compute_event_displace.html - detect event on atom displacement
|
||||
"fep"_compute_fep.html -
|
||||
"force/tally"_compute_tally.html -
|
||||
"fragment/atom"_compute_cluster_atom.html - fragment ID for each atom
|
||||
"global/atom"_compute_global_atom.html -
|
||||
"group/group"_compute_group_group.html - energy/force between two groups of atoms
|
||||
"gyration"_compute_gyration.html - radius of gyration of group of atoms
|
||||
"gyration/chunk"_compute_gyration_chunk.html - radius of gyration for each chunk
|
||||
"heat/flux"_compute_heat_flux.html - heat flux through a group of atoms
|
||||
"heat/flux/tally"_compute_tally.html -
|
||||
"hexorder/atom"_compute_hexorder_atom.html - bond orientational order parameter q6
|
||||
"improper"_compute_improper.html -
|
||||
"hma"_compute_hma.html - harmonically mapped averaging for atomic crystals
|
||||
"improper/local"_compute_improper_local.html - angle of each improper
|
||||
"inertia/chunk"_compute_inertia_chunk.html - inertia tensor for each chunk
|
||||
"ke"_compute_ke.html - translational kinetic energy
|
||||
"ke/atom"_compute_ke_atom.html - kinetic energy for each atom
|
||||
"ke/atom/eff"_compute_ke_atom_eff.html -
|
||||
"ke/eff"_compute_ke_eff.html -
|
||||
"ke/rigid"_compute_ke_rigid.html - translational kinetic energy of rigid bodies
|
||||
"meso/e/atom"_compute_meso_e_atom.html -
|
||||
"meso/rho/atom"_compute_meso_rho_atom.html -
|
||||
"meso/t/atom"_compute_meso_t_atom.html -
|
||||
"msd"_compute_msd.html - mean-squared displacement of group of atoms
|
||||
"msd/chunk"_compute_msd_chunk.html - mean-squared displacement for each chunk
|
||||
"msd/nongauss"_compute_msd_nongauss.html - MSD and non-Gaussian parameter of group of atoms
|
||||
@ -241,77 +219,36 @@ compute"_Commands_compute.html doc page are followed by one or more of
|
||||
"pair/local"_compute_pair_local.html - distance/energy/force of each pairwise interaction
|
||||
"pe"_compute_pe.html - potential energy
|
||||
"pe/atom"_compute_pe_atom.html - potential energy for each atom
|
||||
"pe/mol/tally"_compute_tally.html -
|
||||
"pe/tally"_compute_tally.html -
|
||||
"plasticity/atom"_compute_plasticity_atom.html - Peridynamic plasticity for each atom
|
||||
"pressure"_compute_pressure.html - total pressure and pressure tensor
|
||||
"pressure/cylinder"_compute_pressure_cylinder.html -
|
||||
"pressure/uef"_compute_pressure_uef.html -
|
||||
"property/atom"_compute_property_atom.html - convert atom attributes to per-atom vectors/arrays
|
||||
"property/chunk"_compute_property_chunk.html - extract various per-chunk attributes
|
||||
"property/local"_compute_property_local.html - convert local attributes to localvectors/arrays
|
||||
"ptm/atom"_compute_ptm_atom.html -
|
||||
"property/chunk"_compute_property_chunk.html - extract various per-chunk attributes
|
||||
"rdf"_compute_rdf.html - radial distribution function g(r) histogram of group of atoms
|
||||
"reduce"_compute_reduce.html - combine per-atom quantities into a single global value
|
||||
"reduce/chunk"_compute_reduce_chunk.html - reduce per-atom quantities within each chunk
|
||||
"reduce/region"_compute_reduce.html - same as compute reduce, within a region
|
||||
"rigid/local"_compute_rigid_local.html - extract rigid body attributes
|
||||
"saed"_compute_saed.html -
|
||||
"slice"_compute_slice.html - extract values from global vector or array
|
||||
"smd/contact/radius"_compute_smd_contact_radius.html -
|
||||
"smd/damage"_compute_smd_damage.html -
|
||||
"smd/hourglass/error"_compute_smd_hourglass_error.html -
|
||||
"smd/internal/energy"_compute_smd_internal_energy.html -
|
||||
"smd/plastic/strain"_compute_smd_plastic_strain.html -
|
||||
"smd/plastic/strain/rate"_compute_smd_plastic_strain_rate.html -
|
||||
"smd/rho"_compute_smd_rho.html -
|
||||
"smd/tlsph/defgrad"_compute_smd_tlsph_defgrad.html -
|
||||
"smd/tlsph/dt"_compute_smd_tlsph_dt.html -
|
||||
"smd/tlsph/num/neighs"_compute_smd_tlsph_num_neighs.html -
|
||||
"smd/tlsph/shape"_compute_smd_tlsph_shape.html -
|
||||
"smd/tlsph/strain"_compute_smd_tlsph_strain.html -
|
||||
"smd/tlsph/strain/rate"_compute_smd_tlsph_strain_rate.html -
|
||||
"smd/tlsph/stress"_compute_smd_tlsph_stress.html -
|
||||
"smd/triangle/vertices"_compute_smd_triangle_vertices.html -
|
||||
"smd/triangle/vertices"_compute_smd_triangle_vertices.html -
|
||||
"smd/ulsph/num/neighs"_compute_smd_ulsph_num_neighs.html -
|
||||
"smd/ulsph/strain"_compute_smd_ulsph_strain.html -
|
||||
"smd/ulsph/strain/rate"_compute_smd_ulsph_strain_rate.html -
|
||||
"smd/ulsph/stress"_compute_smd_ulsph_stress.html -
|
||||
"smd/vol"_compute_smd_vol.html -
|
||||
"sna/atom"_compute_sna_atom.html - calculate bispectrum coefficients for each atom
|
||||
"snad/atom"_compute_sna_atom.html - derivative of bispectrum coefficients for each atom
|
||||
"snav/atom"_compute_sna_atom.html - virial contribution from bispectrum coefficients for each atom
|
||||
"spin"_compute_spin.html -
|
||||
"stress/atom"_compute_stress_atom.html - stress tensor for each atom
|
||||
"stress/mop"_compute_stress_mop.html -
|
||||
"stress/mop/profile"_compute_stress_mop.html -
|
||||
"stress/tally"_compute_tally.html -
|
||||
"tdpd/cc/atom"_compute_tdpd_cc_atom.html -
|
||||
"temp"_compute_temp.html - temperature of group of atoms
|
||||
"temp/asphere"_compute_temp_asphere.html - temperature of aspherical particles
|
||||
"temp/body"_compute_temp_body.html - temperature of body particles
|
||||
"temp/chunk"_compute_temp_chunk.html - temperature of each chunk
|
||||
"temp/com"_compute_temp_com.html - temperature after subtracting center-of-mass velocity
|
||||
"temp/cs"_compute_temp_cs.html -
|
||||
"temp/deform"_compute_temp_deform.html - temperature excluding box deformation velocity
|
||||
"temp/deform/eff"_compute_temp_deform_eff.html -
|
||||
"temp/drude"_compute_temp_drude.html -
|
||||
"temp/eff"_compute_temp_eff.html -
|
||||
"temp/partial"_compute_temp_partial.html - temperature excluding one or more dimensions of velocity
|
||||
"temp/profile"_compute_temp_profile.html - temperature excluding a binned velocity profile
|
||||
"temp/ramp"_compute_temp_ramp.html - temperature excluding ramped velocity component
|
||||
"temp/region"_compute_temp_region.html - temperature of a region of atoms
|
||||
"temp/region/eff"_compute_temp_region_eff.html -
|
||||
"temp/rotate"_compute_temp_rotate.html -
|
||||
"temp/sphere"_compute_temp_sphere.html - temperature of spherical particles
|
||||
"temp/uef"_compute_temp_uef.html -
|
||||
"ti"_compute_ti.html - thermodynamic integration free energy values
|
||||
"torque/chunk"_compute_torque_chunk.html - torque applied on each chunk
|
||||
"vacf"_compute_vacf.html - velocity-autocorrelation function of group of atoms
|
||||
"vcm/chunk"_compute_vcm_chunk.html - velocity of center-of-mass for each chunk
|
||||
"voronoi/atom"_compute_voronoi_atom.html - Voronoi volume and neighbors for each atom
|
||||
"xrd"_compute_xrd.html - :ul
|
||||
"voronoi/atom"_compute_voronoi_atom.html - Voronoi volume and neighbors for each atom :ul
|
||||
|
||||
[Restrictions:] none
|
||||
|
||||
|
||||
147
doc/src/compute_hma.txt
Normal file
147
doc/src/compute_hma.txt
Normal file
@ -0,0 +1,147 @@
|
||||
"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c
|
||||
|
||||
:link(lws,http://lammps.sandia.gov)
|
||||
:link(ld,Manual.html)
|
||||
:link(lc,Section_commands.html#comm)
|
||||
|
||||
:line
|
||||
|
||||
compute hma command :h3
|
||||
|
||||
[Syntax:]
|
||||
|
||||
compute ID group-ID hma temp-ID keyword ... :pre
|
||||
|
||||
ID, group-ID are documented in "compute"_compute.html command :l
|
||||
hma = style name of this compute command :l
|
||||
temp-ID = ID of fix that specifies the set temperature during canonical simulation :l
|
||||
keyword = {anharmonic} {u} {p Pharm} {cv} :l
|
||||
{anharmonic} = compute will return anharmonic property values
|
||||
{u} = compute will return potential energy
|
||||
{p} = compute will return pressure. the following keyword must be the difference between the harmonic pressure and lattice pressure as described below
|
||||
{cv} = compute will return the heat capacity :pre
|
||||
:ule
|
||||
|
||||
[Examples:]
|
||||
|
||||
compute 2 all hma 1 u
|
||||
compute 2 all hma 1 anharmonic u p 0.9
|
||||
compute 2 all hma 1 u cv :pre
|
||||
|
||||
|
||||
|
||||
[Description:]
|
||||
|
||||
Define a computation that calculates the properties of a solid (potential
|
||||
energy, pressure or heat capacity), using the harmonically-mapped averaging
|
||||
(HMA) method.
|
||||
This command yields much higher precision than the equivalent compute commands
|
||||
("compute pe"_compute_pe.html, "compute pressure"_compute_pressure.html, etc.)
|
||||
commands during a canonical simulation of an atomic crystal.
|
||||
|
||||
|
||||
In this method, the analytically known harmonic behavior of a crystal is removed from the traditional ensemble
|
||||
averages, which leads to an accurate and precise measurement of the anharmonic contributions without contamination
|
||||
by noise produced by the already-known harmonic behavior.
|
||||
A detailed description of this method can be found in ("Moustafa"_#hma-Moustafa). The potential energy is computed by the formula:
|
||||
|
||||
\begin\{equation\}
|
||||
\left< U\right>_\{HMA\} = \frac\{d(N-1)\}\{2\beta\} + \left< U + \frac\{1\}\{2\} F\bullet\Delta r \right>
|
||||
\end\{equation\}
|
||||
|
||||
where \(N\) is the number of atoms in the system, \(\beta\) is the reciprocal of the thermodynamic temperature, \(d\) is the
|
||||
dimensionality of the system (2 or 3 for 2d/3d), \(F\bullet\Delta r\) is the sum of dot products of the
|
||||
atomic force vectors and displacement (from lattice sites) vectors, and \(U\) is the sum of
|
||||
pair, bond, angle, dihedral, improper, kspace (long-range), and fix energies.
|
||||
|
||||
The pressure is computed by the formula:
|
||||
|
||||
\begin\{equation\}
|
||||
\left< P\right>_\{HMA\} = \Delta \hat P + \left< P_\{vir\} + \frac\{\beta \Delta \hat P - \rho\}\{d(N-1)\} F\bullet\Delta r \right>
|
||||
\end\{equation\}
|
||||
|
||||
where \(\rho\) is the number density of the system, \(\Delta \hat P\) is the
|
||||
difference between theh harmonic and lattice pressure, and \(P_\{vir\}\) is
|
||||
the virial pressure computed as the sum of pair, bond, angle, dihedral,
|
||||
improper, kspace (long-range), and fix contributions to the force on each
|
||||
atom. Although the method will work for any value of \(\Delta \hat P\)
|
||||
specified (use pressure "units"_units.html), the precision of the resultant
|
||||
pressure is sensitive to \(\Delta \hat P\); the precision tends to be
|
||||
best when \(\Delta \hat P\) is the actual the difference between the lattice
|
||||
pressure and harmonic pressure.
|
||||
|
||||
\begin\{equation\}
|
||||
\left<C_V \right>_\{HMA\} = \frac\{d k_B (N-1)\}\{2\} + \beta \left( \left<
|
||||
U_\{HMA\}^2 \right> - \left<U_\{HMA\}\right>^2 \right)/T + \frac\{1\}\{4 T\}
|
||||
\left< F\bullet\Delta r + \Delta r \bullet \Phi \bullet \Delta r \right>
|
||||
\end\{equation\}
|
||||
|
||||
where \(\Phi\) is the Hessian of second derivatives. The compute hma command
|
||||
computes the full expression for \(C_V\) except for the
|
||||
\(\left<U_\{HMA\}^2\right>^2\) in the variance term, which can be obtained by
|
||||
passing the {u} keyword; you must add this extra contribution to the \(C_V\)
|
||||
value reported by this compute. The variance term can cause significant
|
||||
roundoff error when computing \(C_V\). To address this, the {anharmonic}
|
||||
keyword can be passed and/or the output format can be speicified with more
|
||||
digits.
|
||||
|
||||
thermo_modify format float '%22.15e' :pre
|
||||
|
||||
The {anharmonic} keyword will instruct the compute to return anharmonic
|
||||
properties rather than the full properties (lattice, harmonic and anharmonic).
|
||||
When using this keyword, the compute must be first active (it must be included
|
||||
via a "thermo_style custom"_thermo_style.html command) while the atoms are
|
||||
still at their lattice sites (before equilibration).
|
||||
|
||||
The temp-ID specified with compute hma command should be same as the fix-ID of Nose-Hoover ("fix nvt"_fix_nh.html) or
|
||||
Berendsen ("fix temp/berendsen"_fix_temp_berendsen.html) thermostat used for the simulation. While using this command, Langevin thermostat
|
||||
("fix langevin"_fix_langevin.html)
|
||||
should be avoided as its extra forces interfere with the HMA implementation.
|
||||
|
||||
|
||||
|
||||
NOTE: Compute hma command should be used right after the energy minimization, when the atoms are at their lattice sites.
|
||||
The simulation should not be started before this command has been used in the input script.
|
||||
|
||||
|
||||
The following example illustrates the placement of this command in the input script:
|
||||
|
||||
|
||||
min_style cg
|
||||
minimize 1e-35 1e-15 50000 500000
|
||||
compute 1 all hma thermostatid u
|
||||
fix thermostatid all nvt temp 600.0 600.0 100.0 :pre
|
||||
|
||||
|
||||
|
||||
NOTE: Compute hma should be used when the atoms of the solid do not diffuse. Diffusion will reduce the precision in the potential energy computation.
|
||||
|
||||
|
||||
NOTE: The "fix_modify energy yes"_fix_modify.html command must also be specified if a fix is to contribute potential energy to this command.
|
||||
|
||||
:line
|
||||
|
||||
[Output info:]
|
||||
|
||||
This compute calculates a global vector that includes the n properties
|
||||
requested as arguments to the command (the potential energy, pressure or heat
|
||||
capacity). The elements of the vector can be accessed by indices 1-n by any
|
||||
command that uses global vector values as input. See the "Howto
|
||||
output"_Howto_output.html doc page for an overview of LAMMPS output options.
|
||||
|
||||
The vector values calculated by this compute are "extensive". The
|
||||
scalar value will be in energy "units"_units.html.
|
||||
|
||||
[Restrictions:] Usage restricted to canonical (NVT) ensemble simulation only.
|
||||
|
||||
[Related commands:]
|
||||
|
||||
"compute pe"_compute_pe.html, "compute pressure"_compute_pressure.html
|
||||
|
||||
[Default:] none
|
||||
|
||||
:line
|
||||
|
||||
:link(hma-Moustafa)
|
||||
[(Moustafa)] Sabry G. Moustafa, Andrew J. Schultz, and David A. Kofke, {Very fast averaging of thermal properties of crystals by molecular simulation},
|
||||
"Phys. Rev. E \[92\], 043303 (2015)"_https://link.aps.org/doi/10.1103/PhysRevE.92.043303
|
||||
@ -6,7 +6,6 @@ Computes :h1
|
||||
:maxdepth: 1
|
||||
|
||||
compute_ackland_atom
|
||||
compute_adf
|
||||
compute_angle
|
||||
compute_angle_local
|
||||
compute_angmom_chunk
|
||||
@ -16,7 +15,6 @@ Computes :h1
|
||||
compute_bond_local
|
||||
compute_centro_atom
|
||||
compute_chunk_atom
|
||||
compute_chunk_spread_atom
|
||||
compute_cluster_atom
|
||||
compute_cna_atom
|
||||
compute_cnp_atom
|
||||
@ -46,6 +44,7 @@ Computes :h1
|
||||
compute_gyration_chunk
|
||||
compute_heat_flux
|
||||
compute_hexorder_atom
|
||||
compute_hma
|
||||
compute_improper
|
||||
compute_improper_local
|
||||
compute_inertia_chunk
|
||||
@ -68,15 +67,12 @@ Computes :h1
|
||||
compute_pe_atom
|
||||
compute_plasticity_atom
|
||||
compute_pressure
|
||||
compute_pressure_cylinder
|
||||
compute_pressure_uef
|
||||
compute_property_atom
|
||||
compute_property_chunk
|
||||
compute_property_local
|
||||
compute_ptm_atom
|
||||
compute_rdf
|
||||
compute_reduce
|
||||
compute_reduce_chunk
|
||||
compute_rigid_local
|
||||
compute_saed
|
||||
compute_slice
|
||||
@ -94,7 +90,7 @@ Computes :h1
|
||||
compute_smd_tlsph_strain
|
||||
compute_smd_tlsph_strain_rate
|
||||
compute_smd_tlsph_stress
|
||||
compute_smd_triangle_vertices
|
||||
compute_smd_triangle_mesh_vertices
|
||||
compute_smd_ulsph_num_neighs
|
||||
compute_smd_ulsph_strain
|
||||
compute_smd_ulsph_strain_rate
|
||||
@ -103,7 +99,6 @@ Computes :h1
|
||||
compute_sna_atom
|
||||
compute_spin
|
||||
compute_stress_atom
|
||||
compute_stress_mop
|
||||
compute_tally
|
||||
compute_tdpd_cc_atom
|
||||
compute_temp
|
||||
|
||||
470
src/compute_hma.cpp
Normal file
470
src/compute_hma.cpp
Normal file
@ -0,0 +1,470 @@
|
||||
/* ----------------------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
http://lammps.sandia.gov, Sandia National Laboratories
|
||||
Steve Plimpton, sjplimp@sandia.gov
|
||||
|
||||
Copyright (2003) Sandia Corporation. Under the terms of Contract
|
||||
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
|
||||
certain rights in this software. This software is distributed under
|
||||
the GNU General Public License.
|
||||
|
||||
See the README file in the top-level LAMMPS directory.
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#include <mpi.h>
|
||||
#include <math.h>
|
||||
#include <string.h>
|
||||
#include "compute_hma.h"
|
||||
#include "atom.h"
|
||||
#include "update.h"
|
||||
#include "force.h"
|
||||
#include "pair.h"
|
||||
#include "bond.h"
|
||||
#include "angle.h"
|
||||
#include "dihedral.h"
|
||||
#include "improper.h"
|
||||
#include "kspace.h"
|
||||
#include "group.h"
|
||||
#include "domain.h"
|
||||
#include "modify.h"
|
||||
#include "fix.h"
|
||||
#include "fix_store.h"
|
||||
#include "memory.h"
|
||||
#include "error.h"
|
||||
#include "comm.h"
|
||||
#include "neighbor.h"
|
||||
#include "neigh_request.h"
|
||||
#include "neigh_list.h"
|
||||
|
||||
#include <vector>
|
||||
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
ComputeHMA::ComputeHMA(LAMMPS *lmp, int narg, char **arg) :
|
||||
Compute(lmp, narg, arg), id_temp(NULL), deltaR(NULL)
|
||||
{
|
||||
if (narg < 4) error->all(FLERR,"Illegal compute hma command");
|
||||
if (igroup) error->all(FLERR,"Compute hma must use group all");
|
||||
if (strcmp(arg[3],"NULL") == 0) {error->all(FLERR,"fix ID specifying the set temperature of canonical simulation is required");}
|
||||
else {
|
||||
int n = strlen(arg[3]) + 1;
|
||||
id_temp = new char[n];
|
||||
strcpy(id_temp,arg[3]);
|
||||
}
|
||||
|
||||
create_attribute = 1;
|
||||
extscalar = 1;
|
||||
timeflag = 1;
|
||||
|
||||
// (from compute displace/atom) create a new fix STORE style
|
||||
// our new fix's id (id_fix)= compute-ID + COMPUTE_STORE
|
||||
// our new fix's group = same as compute group
|
||||
|
||||
int n = strlen(id) + strlen("_COMPUTE_STORE") + 1;
|
||||
id_fix = new char[n];
|
||||
strcpy(id_fix,id);
|
||||
strcat(id_fix,"_COMPUTE_STORE");
|
||||
|
||||
char **newarg = new char*[6];
|
||||
newarg[0] = id_fix;
|
||||
newarg[1] = group->names[igroup];
|
||||
newarg[2] = (char *) "STORE";
|
||||
newarg[3] = (char *) "peratom";
|
||||
newarg[4] = (char *) "1";
|
||||
newarg[5] = (char *) "3";
|
||||
modify->add_fix(6,newarg);
|
||||
fix = (FixStore *) modify->fix[modify->nfix-1];
|
||||
|
||||
delete [] newarg;
|
||||
|
||||
// calculate xu,yu,zu for fix store array
|
||||
// skip if reset from restart file
|
||||
|
||||
if (fix->restart_reset) fix->restart_reset = 0;
|
||||
else {
|
||||
double **xoriginal = fix->astore;
|
||||
double **x = atom->x;
|
||||
imageint *image = atom->image;
|
||||
int nlocal = atom->nlocal;
|
||||
|
||||
for (int i = 0; i < nlocal; i++)
|
||||
domain->unmap(x[i],image[i],xoriginal[i]);
|
||||
}
|
||||
|
||||
vector_flag = 1;
|
||||
extvector = -1;
|
||||
comm_forward = 0; // 3 if 2nd derivative needed
|
||||
|
||||
computeU = computeP = computeCv = -1;
|
||||
returnAnharmonic = 0;
|
||||
std::vector<int> extvec;
|
||||
for (int iarg=4; iarg<narg; iarg++) {
|
||||
if (!strcasecmp(arg[iarg], "u")) {
|
||||
computeU = extvec.size();
|
||||
extvec.push_back(1);
|
||||
}
|
||||
else if (!strcasecmp(arg[iarg], "p")) {
|
||||
if (iarg+2 > narg) error->all(FLERR,"Illegal fix hma command");
|
||||
computeP = extvec.size();
|
||||
deltaPcap = force->numeric(FLERR, arg[iarg+1]);
|
||||
extvec.push_back(0);
|
||||
iarg++;
|
||||
}
|
||||
else if (!strcasecmp(arg[iarg], "cv")) {
|
||||
computeCv = extvec.size();
|
||||
comm_forward = 3;
|
||||
extvec.push_back(1);
|
||||
}
|
||||
else if (!strcasecmp(arg[iarg], "anharmonic")) {
|
||||
returnAnharmonic = -1;
|
||||
}
|
||||
else {
|
||||
error->all(FLERR,"Illegal fix hma command");
|
||||
}
|
||||
}
|
||||
|
||||
if (extvec.size() == 0) {
|
||||
error->all(FLERR,"Illegal fix hma command");
|
||||
}
|
||||
size_vector = extvec.size();
|
||||
memory->create(vector, size_vector, "hma::vector");
|
||||
extlist = new int[size_vector];
|
||||
for (int i=0; i<size_vector; i++) {
|
||||
extlist[i] = extvec[i];
|
||||
}
|
||||
|
||||
if (computeU>-1 || computeCv>-1) {
|
||||
peflag = 1;
|
||||
}
|
||||
if (computeP>-1) {
|
||||
pressflag = 1;
|
||||
}
|
||||
|
||||
nmax = 0;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
ComputeHMA::~ComputeHMA()
|
||||
{
|
||||
// check nfix in case all fixes have already been deleted
|
||||
if (modify->nfix) modify->delete_fix(id_fix);
|
||||
|
||||
delete [] id_fix;
|
||||
delete [] id_temp;
|
||||
delete [] extlist;
|
||||
memory->destroy(vector);
|
||||
memory->destroy(deltaR);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void ComputeHMA::init() {
|
||||
int irequest = neighbor->request(this,instance_me);
|
||||
neighbor->requests[irequest]->pair = 0;
|
||||
neighbor->requests[irequest]->compute = 1;
|
||||
neighbor->requests[irequest]->occasional = 1;
|
||||
}
|
||||
|
||||
void ComputeHMA::init_list(int id, NeighList *ptr) {
|
||||
list = ptr;
|
||||
}
|
||||
|
||||
void ComputeHMA::setup()
|
||||
{
|
||||
int dummy=0;
|
||||
int ifix = modify->find_fix(id_temp);
|
||||
if (ifix < 0) error->all(FLERR,"Could not find compute hma temperature ID");
|
||||
double * temperat = (double *) modify->fix[ifix]->extract("t_target",dummy);
|
||||
if (temperat==NULL) error->all(FLERR,"Could not find compute hma temperature ID");
|
||||
finaltemp = * temperat;
|
||||
|
||||
|
||||
// set fix which stores original atom coords
|
||||
|
||||
int ifix2 = modify->find_fix(id_fix);
|
||||
if (ifix2 < 0) error->all(FLERR,"Could not find compute hma ID");
|
||||
fix = (FixStore *) modify->fix[ifix2];
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void ComputeHMA::compute_vector()
|
||||
{
|
||||
invoked_vector = update->ntimestep;
|
||||
|
||||
// dx,dy,dz = displacement of atom from original position
|
||||
// original unwrapped position is stored by fix
|
||||
// for triclinic, need to unwrap current atom coord via h matrix
|
||||
|
||||
// grow deltaR array if necessary
|
||||
|
||||
if (comm_forward>0 && atom->nmax > nmax) {
|
||||
memory->destroy(deltaR);
|
||||
nmax = atom->nmax;
|
||||
memory->create(deltaR,nmax,3,"hma:deltaR");
|
||||
}
|
||||
|
||||
double **xoriginal = fix->astore;
|
||||
double fdr = 0.0;
|
||||
double **x = atom->x;
|
||||
double **f = atom->f;
|
||||
|
||||
imageint *image = atom->image;
|
||||
int nlocal = atom->nlocal;
|
||||
|
||||
double *h = domain->h;
|
||||
double xprd = domain->xprd;
|
||||
double yprd = domain->yprd;
|
||||
double zprd = domain->zprd;
|
||||
|
||||
double u = 0.0;
|
||||
if (computeU>-1 || computeCv>-1) {
|
||||
if (force->pair) u += force->pair->eng_vdwl + force->pair->eng_coul;
|
||||
if (force->bond) u += force->bond->energy;
|
||||
if (force->angle) u += force->angle->energy;
|
||||
if (force->dihedral) u += force->dihedral->energy;
|
||||
if (force->improper) u += force->improper->energy;
|
||||
}
|
||||
|
||||
int dimension = domain->dimension;
|
||||
double p = 0, vol = 0;
|
||||
if (computeP>-1) {
|
||||
p = virial_compute(dimension);
|
||||
vol = xprd * yprd;
|
||||
if (dimension == 3) vol *= zprd;
|
||||
p *= force->nktv2p / (dimension*vol);
|
||||
if (returnAnharmonic == -1) {
|
||||
pLat = p;
|
||||
}
|
||||
}
|
||||
|
||||
if (domain->triclinic == 0) {
|
||||
for (int i = 0; i < nlocal; i++) {
|
||||
int xbox = (image[i] & IMGMASK) - IMGMAX;
|
||||
int ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
|
||||
int zbox = (image[i] >> IMG2BITS) - IMGMAX;
|
||||
double dx = x[i][0] + xbox*xprd - xoriginal[i][0];
|
||||
double dy = x[i][1] + ybox*yprd - xoriginal[i][1];
|
||||
double dz = x[i][2] + zbox*zprd - xoriginal[i][2];
|
||||
if (comm_forward>0) {
|
||||
deltaR[i][0] = dx;
|
||||
deltaR[i][1] = dy;
|
||||
deltaR[i][2] = dz;
|
||||
}
|
||||
fdr += dx*f[i][0] + dy*f[i][1] + dz*f[i][2];
|
||||
}
|
||||
}
|
||||
else {
|
||||
for (int i = 0; i < nlocal; i++) {
|
||||
int xbox = (image[i] & IMGMASK) - IMGMAX;
|
||||
int ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
|
||||
int zbox = (image[i] >> IMG2BITS) - IMGMAX;
|
||||
double dx = x[i][0] + h[0]*xbox + h[5]*ybox + h[4]*zbox - xoriginal[i][0];
|
||||
double dy = x[i][1] + h[1]*ybox + h[3]*zbox - xoriginal[i][1];
|
||||
double dz = x[i][2] + h[2]*zbox - xoriginal[i][2];
|
||||
if (comm_forward>0) {
|
||||
deltaR[i][0] = dx;
|
||||
deltaR[i][1] = dy;
|
||||
deltaR[i][2] = dz;
|
||||
}
|
||||
fdr += ((dx*f[i][0])+(dy*f[i][1])+(dz*f[i][2]));
|
||||
}
|
||||
}
|
||||
|
||||
double phiSum = 0.0;
|
||||
if (computeCv>-1) {
|
||||
comm->forward_comm_compute(this);
|
||||
int *type = atom->type;
|
||||
double** cutsq = force->pair->cutsq;
|
||||
if (force->pair) {
|
||||
double **x = atom->x;
|
||||
double **f = atom->f;
|
||||
int *type = atom->type;
|
||||
int nlocal = atom->nlocal;
|
||||
double *special_lj = force->special_lj;
|
||||
double *special_coul = force->special_coul;
|
||||
int newton_pair = force->newton_pair;
|
||||
|
||||
if (update->firststep == update->ntimestep) neighbor->build_one(list,1);
|
||||
else neighbor->build_one(list);
|
||||
int inum = list->inum;
|
||||
int *ilist = list->ilist;
|
||||
int *numneigh = list->numneigh;
|
||||
int **firstneigh = list->firstneigh;
|
||||
|
||||
for (int ii = 0; ii < inum; ii++) {
|
||||
int i = ilist[ii];
|
||||
double fac = (newton_pair || i < nlocal) ? 1.0 : 0.5;
|
||||
double* ix = x[i];
|
||||
int itype = type[i];
|
||||
int *jlist = firstneigh[i];
|
||||
int jnum = numneigh[i];
|
||||
double *idr = deltaR[i];
|
||||
for (int jj = 0; jj < jnum; jj++) {
|
||||
int j = jlist[jj];
|
||||
if (!newton_pair && j>=nlocal) fac -= 0.5;
|
||||
double factor_lj = special_lj[sbmask(j)];
|
||||
double factor_coul = special_coul[sbmask(j)];
|
||||
j &= NEIGHMASK;
|
||||
double* jx = x[j];
|
||||
double delr[3];
|
||||
delr[0] = ix[0] - jx[0];
|
||||
delr[1] = ix[1] - jx[1];
|
||||
delr[2] = ix[2] - jx[2];
|
||||
double rsq = delr[0]*delr[0] + delr[1]*delr[1] + delr[2]*delr[2];
|
||||
int jtype = type[j];
|
||||
if (rsq < cutsq[itype][jtype]) {
|
||||
double* jdr = deltaR[j];
|
||||
double fforce, d2u[6];
|
||||
force->pair->single2(i, j, itype, jtype, rsq, delr, factor_coul, factor_lj, fforce, d2u);
|
||||
int m = 0;
|
||||
for (int k=0; k<3; k++) {
|
||||
double a = fac;
|
||||
for (int l=k; l<3; l++) {
|
||||
phiSum += a*(idr[k]*jdr[l]+jdr[k]*idr[l])*d2u[m];
|
||||
phiSum -= a*(idr[k]*idr[l]*d2u[m] + jdr[k]*jdr[l]*d2u[m]);
|
||||
m++;
|
||||
if (k==l) a *= 2;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
//adding PE for this processor
|
||||
|
||||
//adding the energies of all processors
|
||||
|
||||
double fdrTotal;
|
||||
MPI_Allreduce(&fdr,&fdrTotal,1,MPI_DOUBLE,MPI_SUM,world);
|
||||
double uTotal;
|
||||
if (computeU>-1 || computeCv>-1) {
|
||||
MPI_Allreduce(&u,&uTotal,1,MPI_DOUBLE,MPI_SUM,world);
|
||||
if (returnAnharmonic == -1) {
|
||||
uLat = uTotal;
|
||||
}
|
||||
if (computeU>-1) {
|
||||
if (returnAnharmonic) {
|
||||
vector[computeU] = uTotal - uLat + 0.5*fdrTotal;
|
||||
}
|
||||
else {
|
||||
vector[computeU] = uTotal + 0.5*fdrTotal + 0.5*dimension*(atom->natoms - 1)*force->boltz*finaltemp;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
if (computeP>-1) {
|
||||
double fv = ((deltaPcap)-(force->boltz*finaltemp*force->nktv2p*atom->natoms/vol))/(force->boltz*finaltemp*dimension*(atom->natoms - 1));
|
||||
if (returnAnharmonic) {
|
||||
vector[computeP] = p - pLat + (fv*fdrTotal);
|
||||
}
|
||||
else {
|
||||
vector[computeP] = p + (fv*fdrTotal) + deltaPcap;
|
||||
}
|
||||
}
|
||||
|
||||
if (computeCv>-1) {
|
||||
if (computeU==-1) MPI_Allreduce(&u,&uTotal,1,MPI_DOUBLE,MPI_SUM,world);
|
||||
double buTot;
|
||||
if (returnAnharmonic) {
|
||||
buTot = (uTotal - uLat + 0.5*fdrTotal)/finaltemp;
|
||||
}
|
||||
else {
|
||||
buTot = (uTotal + 0.5*fdrTotal)/finaltemp + 0.5*dimension*(atom->natoms - 1)*force->boltz;
|
||||
}
|
||||
double one = -0.25*(fdr + phiSum)/finaltemp;
|
||||
double Cv;
|
||||
MPI_Allreduce(&one,&Cv,1,MPI_DOUBLE,MPI_SUM,world);
|
||||
vector[computeCv] = Cv + buTot*buTot;
|
||||
if (!returnAnharmonic) {
|
||||
vector[computeCv] += 0.5*dimension*(atom->natoms-1);
|
||||
}
|
||||
}
|
||||
if (returnAnharmonic == -1) {
|
||||
returnAnharmonic = 1;
|
||||
}
|
||||
}
|
||||
|
||||
double ComputeHMA::virial_compute(int n) {
|
||||
double v = 0;
|
||||
|
||||
// sum contributions to virial from forces and fixes
|
||||
|
||||
if (force->pair) v += sumVirial(n, force->pair->virial);
|
||||
if (force->bond) v += sumVirial(n, force->bond->virial);
|
||||
if (force->angle) v += sumVirial(n, force->angle->virial);
|
||||
if (force->dihedral) v += sumVirial(n, force->dihedral->virial);
|
||||
if (force->improper) v += sumVirial(n, force->improper->virial);
|
||||
for (int i = 0; i < modify->nfix; i++)
|
||||
if (modify->fix[i]->thermo_virial) v += sumVirial(n, modify->fix[i]->virial);
|
||||
|
||||
// sum virial across procs
|
||||
|
||||
double virial;
|
||||
MPI_Allreduce(&v,&virial,1,MPI_DOUBLE,MPI_SUM,world);
|
||||
|
||||
// KSpace virial contribution is already summed across procs
|
||||
|
||||
if (force->kspace)
|
||||
for (int i = 0; i < n; i++) virial += force->kspace->virial[i];
|
||||
return virial;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
int ComputeHMA::pack_forward_comm(int n, int *list, double *buf,
|
||||
int pbc_flag, int *pbc)
|
||||
{
|
||||
double **xoriginal = fix->astore;
|
||||
imageint *image = atom->image;
|
||||
double **x = atom->x;
|
||||
double *h = domain->h;
|
||||
double xprd = domain->xprd;
|
||||
double yprd = domain->yprd;
|
||||
double zprd = domain->zprd;
|
||||
|
||||
int m = 0;
|
||||
for (int ii = 0; ii < n; ii++) {
|
||||
int i = list[ii];
|
||||
buf[m++] = deltaR[i][0];
|
||||
buf[m++] = deltaR[i][1];
|
||||
buf[m++] = deltaR[i][2];
|
||||
}
|
||||
return m;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void ComputeHMA::unpack_forward_comm(int n, int first, double *buf)
|
||||
{
|
||||
double **xoriginal = fix->astore;
|
||||
int i,m,last;
|
||||
|
||||
m = 0;
|
||||
last = first + n;
|
||||
for (i = first; i < last; i++) {
|
||||
deltaR[i][0] = buf[m++];
|
||||
deltaR[i][1] = buf[m++];
|
||||
deltaR[i][2] = buf[m++];
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
initialize one atom's storage values, called when atom is created
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void ComputeHMA::set_arrays(int i)
|
||||
{
|
||||
double **xoriginal = fix->astore;
|
||||
double **x = atom->x;
|
||||
xoriginal[i][0] = x[i][0];
|
||||
xoriginal[i][1] = x[i][1];
|
||||
xoriginal[i][2] = x[i][2];
|
||||
}
|
||||
66
src/compute_hma.h
Normal file
66
src/compute_hma.h
Normal file
@ -0,0 +1,66 @@
|
||||
/* -*- c++ -*- ----------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
http://lammps.sandia.gov, Sandia National Laboratories
|
||||
Steve Plimpton, sjplimp@sandia.gov
|
||||
|
||||
Copyright (2003) Sandia Corporation. Under the terms of Contract
|
||||
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
|
||||
certain rights in this software. This software is distributed under
|
||||
the GNU General Public License.
|
||||
|
||||
See the README file in the top-level LAMMPS directory.
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#ifdef COMPUTE_CLASS
|
||||
|
||||
ComputeStyle(HMA,ComputeHMA)
|
||||
|
||||
#else
|
||||
|
||||
#ifndef LMP_COMPUTE_HMA_H
|
||||
#define LMP_COMPUTE_HMA_H
|
||||
|
||||
#include "compute.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class ComputeHMA : public Compute {
|
||||
public:
|
||||
ComputeHMA(class LAMMPS *, int, char **);
|
||||
~ComputeHMA();
|
||||
void setup();
|
||||
void init();
|
||||
void init_list(int, class NeighList *);
|
||||
void compute_vector();
|
||||
void set_arrays(int);
|
||||
int pack_forward_comm(int, int *, double *, int, int *);
|
||||
void unpack_forward_comm(int, int, double *);
|
||||
|
||||
private:
|
||||
int nmax;
|
||||
int atomsingroup;
|
||||
char *id_fix;
|
||||
char *id_temp;
|
||||
double finaltemp;
|
||||
class FixStore *fix;
|
||||
double boltz, nktv2p, inv_volume;
|
||||
double deltaPcap;
|
||||
double virial_compute(int);
|
||||
static double sumVirial(int n, double* v) {
|
||||
double x = 0;
|
||||
for (int i=0; i<n; i++) x += v[i];
|
||||
return x;
|
||||
}
|
||||
int computeU, computeP, computeCv;
|
||||
class NeighList *list; // half neighbor list
|
||||
double **deltaR;
|
||||
int returnAnharmonic;
|
||||
double uLat, pLat;
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
#endif
|
||||
#endif
|
||||
|
||||
|
||||
Reference in New Issue
Block a user