diff --git a/doc/src/compute.txt b/doc/src/compute.txt index 857795ffe5..676c00d7a5 100644 --- a/doc/src/compute.txt +++ b/doc/src/compute.txt @@ -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 diff --git a/doc/src/compute_hma.txt b/doc/src/compute_hma.txt new file mode 100644 index 0000000000..f71ddc9e2e --- /dev/null +++ b/doc/src/compute_hma.txt @@ -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_\{HMA\} = \frac\{d k_B (N-1)\}\{2\} + \beta \left( \left< +U_\{HMA\}^2 \right> - \left^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^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 diff --git a/doc/src/computes.txt b/doc/src/computes.txt index 926b8da222..c636633b38 100644 --- a/doc/src/computes.txt +++ b/doc/src/computes.txt @@ -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 diff --git a/src/compute_hma.cpp b/src/compute_hma.cpp new file mode 100644 index 0000000000..ac976080ce --- /dev/null +++ b/src/compute_hma.cpp @@ -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 +#include +#include +#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 + + +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 extvec; + for (int iarg=4; iarg 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-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]; +} diff --git a/src/compute_hma.h b/src/compute_hma.h new file mode 100644 index 0000000000..f40103c3c4 --- /dev/null +++ b/src/compute_hma.h @@ -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