tx258/lammps-gran-kokkos
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge pull request #2233 from charlessievers/OptimizedDynamicalMatrix

Browse Source 
Dynamical_matrix and third_order features+update
This commit is contained in:
Axel Kohlmeyer
2022-01-31 12:57:28 -05:00
committed by GitHub
parent a1d186b5fa f537a31c19
commit a124ee4f44
15 changed files with 1593 additions and 425 deletions
Download Patch File Download Diff File Expand all files Collapse all files

6
cmake/Modules/Packages/KOKKOS.cmake
View File

@ -109,6 +109,12 @@ if(PKG_KSPACE)
endif()
endif()
if(PKG_PHONON)
list(APPEND KOKKOS_PKG_SOURCES ${KOKKOS_PKG_SOURCES_DIR}/dynamical_matrix_kokkos.cpp)
list(APPEND KOKKOS_PKG_SOURCES ${KOKKOS_PKG_SOURCES_DIR}/third_order_kokkos.cpp)
endif()
set_property(GLOBAL PROPERTY "KOKKOS_PKG_SOURCES" "${KOKKOS_PKG_SOURCES}")
# detects styles which have KOKKOS version

4
doc/src/Commands_all.rst
View File

@ -47,7 +47,7 @@ An alphabetic list of all general LAMMPS commands.
* :doc:`displace_atoms <displace_atoms>`
* :doc:`dump <dump>`
* :doc:`dump_modify <dump_modify>`
* :doc:`dynamical_matrix <dynamical_matrix>`
* :doc:`dynamical_matrix (k) <dynamical_matrix>`
* :doc:`echo <echo>`
* :doc:`fix <fix>`
* :doc:`fix_modify <fix_modify>`
@ -117,7 +117,7 @@ An alphabetic list of all general LAMMPS commands.
* :doc:`thermo <thermo>`
* :doc:`thermo_modify <thermo_modify>`
* :doc:`thermo_style <thermo_style>`
* :doc:`third_order <third_order>`
* :doc:`third_order (k) <third_order>`
* :doc:`timer <timer>`
* :doc:`timestep <timestep>`
* :doc:`uncompute <uncompute>`

9
doc/src/dynamical_matrix.rst
View File

@ -1,8 +1,11 @@
.. index:: dynamical_matrix
.. index:: dynamical_matrix/kk
dynamical_matrix command
========================
Accelerator Variants: dynamical_matrix/kk
Syntax
""""""
@ -56,6 +59,12 @@ If the style eskm is selected, the dynamical matrix will be in units of
inverse squared femtoseconds. These units will then conveniently leave
frequencies in THz.
----------
.. include:: accel_styles.rst
----------
Restrictions
""""""""""""

9
doc/src/third_order.rst
View File

@ -1,8 +1,11 @@
.. index:: third_order
.. index:: third_order/kk
third_order command
===================
Accelerator Variants: third_order/kk
Syntax
""""""
@ -49,6 +52,12 @@ If the style eskm is selected, the tensor will be using energy units of 10 J/mol
These units conform to eskm style from the dynamical_matrix command, which
will simplify operations using dynamical matrices with third order tensors.
----------
.. include:: accel_styles.rst
----------
Restrictions
""""""""""""

4
src/Depend.sh
View File

@ -135,6 +135,10 @@ if (test $1 = "PYTHON") then
depend ML-IAP
fi
if (test $1 = "PHONON") then
depend KOKKOS
fi
if (test $1 = "RIGID") then
depend KOKKOS
depend OPENMP

4
src/KOKKOS/Install.sh
View File

@ -108,6 +108,8 @@ action dihedral_opls_kokkos.cpp dihedral_opls.cpp
action dihedral_opls_kokkos.h dihedral_opls.h
action domain_kokkos.cpp
action domain_kokkos.h
action dynamical_matrix_kokkos.cpp dynamical_matrix.cpp
action dynamical_matrix_kokkos.h dynamical_matrix.h
action fftdata_kokkos.h fft3d.h
action fft3d_kokkos.cpp fft3d.cpp
action fft3d_kokkos.h fft3d.h
@ -312,6 +314,8 @@ action remap_kokkos.cpp remap.cpp
action remap_kokkos.h remap.h
action sna_kokkos.h sna.h
action sna_kokkos_impl.h sna.cpp
action third_order_kokkos.cpp dynamical_matrix.cpp
action third_order_kokkos.h dynamical_matrix.h
action verlet_kokkos.cpp
action verlet_kokkos.h

353
src/KOKKOS/dynamical_matrix_kokkos.cpp Normal file
View File

@ -0,0 +1,353 @@
// clang-format off
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Charlie Sievers (UC Davis), charliesievers at cox.net
------------------------------------------------------------------------- */
#include "dynamical_matrix_kokkos.h"
#include "angle.h"
#include "atom_kokkos.h"
#include "atom_masks.h"
#include "bond.h"
#include "comm.h"
#include "compute.h"
#include "dihedral.h"
#include "domain.h"
#include "error.h"
#include "finish.h"
#include "force.h"
#include "group.h"
#include "improper.h"
#include "kokkos.h"
#include "kspace.h"
#include "memory.h"
#include "modify.h"
#include "neighbor.h"
#include "pair.h"
#include "timer.h"
#include "update.h"
#include <cmath>
#include <cstring>
#include <algorithm>
using namespace LAMMPS_NS;
enum{REGULAR,ESKM};
template<class ViewA, class ViewB>
struct ForceAdder {
ViewA a;
ViewB b;
ForceAdder(const ViewA& a_, const ViewB& b_):a(a_),b(b_) {}
KOKKOS_INLINE_FUNCTION
void operator() (const int& i) const {
a(i,0) += b(i,0);
a(i,1) += b(i,1);
a(i,2) += b(i,2);
}
};
/* ---------------------------------------------------------------------- */
template<class View>
struct Zero {
View v;
Zero(const View &v_):v(v_) {}
KOKKOS_INLINE_FUNCTION
void operator()(const int &i) const {
v(i,0) = 0;
v(i,1) = 0;
v(i,2) = 0;
}
};
/* ---------------------------------------------------------------------- */
DynamicalMatrixKokkos::DynamicalMatrixKokkos(LAMMPS *lmp) : DynamicalMatrix(lmp)
{
atomKK = (AtomKokkos *) atom;
}
/* ---------------------------------------------------------------------- */
DynamicalMatrixKokkos::~DynamicalMatrixKokkos()
{
}
/* ---------------------------------------------------------------------- */
void DynamicalMatrixKokkos::command(int narg, char **arg)
{
atomKK->sync(Host, X_MASK|RMASS_MASK|TYPE_MASK);
DynamicalMatrix::command(narg, arg);
}
/* ----------------------------------------------------------------------
setup without output or one-time post-init setup
flag = 0 = just force calculation
flag = 1 = reneighbor and force calculation
------------------------------------------------------------------------- */
void DynamicalMatrixKokkos::setup()
{
lmp->kokkos->auto_sync = 1;
// setup domain, communication and neighboring
// acquire ghosts
// build neighbor lists
if (triclinic) domain->x2lamda(atom->nlocal);
domain->pbc();
domain->reset_box();
comm->setup();
if (neighbor->style) neighbor->setup_bins();
comm->exchange();
comm->borders();
if (triclinic) domain->lamda2x(atom->nlocal+atom->nghost);
domain->image_check();
domain->box_too_small_check();
neighbor->build(1);
// compute all forces
eflag=0;
vflag=0;
if (force->kspace) {
force->kspace->setup();
}
update_force();
if (pair_compute_flag) {
atomKK->sync(force->pair->execution_space,force->pair->datamask_read);
force->pair->compute(eflag,vflag);
atomKK->modified(force->pair->execution_space,force->pair->datamask_modify);
}
else if (force->pair) force->pair->compute_dummy(eflag,vflag);
update->setupflag = 0;
lmp->kokkos->auto_sync = 0;
//if all then skip communication groupmap population
if (gcount == atom->natoms)
for (bigint i=0; i<atom->natoms; i++)
groupmap[i] = i;
else
create_groupmap();
}
/* ----------------------------------------------------------------------
evaluate potential energy and forces
may migrate atoms due to reneighboring
return new energy, which should include nextra_global dof
return negative gradient stored in atom->f
return negative gradient for nextra_global dof in fextra
------------------------------------------------------------------------- */
void DynamicalMatrixKokkos::update_force()
{
int n_pre_force = modify->n_pre_force;
int n_pre_reverse = modify->n_pre_reverse;
int n_post_force = modify->n_post_force;
lmp->kokkos->auto_sync = 0;
f_merge_copy = DAT::t_f_array("DynamicalMatrixKokkos::f_merge_copy",atomKK->k_f.extent(0));
atomKK->modified(Host,X_MASK);
atomKK->sync(Device,X_MASK);
force_clear();
neighbor->ago = 0;
if ((modify->get_fix_by_id("package_intel")) ? true : false)
neighbor->decide();
if (n_pre_force) {
modify->pre_force(vflag);
timer->stamp(Timer::MODIFY);
}
bool execute_on_host = false;
unsigned int datamask_read_device = 0;
unsigned int datamask_modify_device = 0;
unsigned int datamask_read_host = 0;
if (pair_compute_flag) {
if (force->pair->execution_space==Host) {
execute_on_host = true;
datamask_read_host |= force->pair->datamask_read;
datamask_modify_device |= force->pair->datamask_modify;
} else {
datamask_read_device |= force->pair->datamask_read;
datamask_modify_device |= force->pair->datamask_modify;
}
}
if (atomKK->molecular && force->bond) {
if (force->bond->execution_space==Host) {
execute_on_host = true;
datamask_read_host |= force->bond->datamask_read;
datamask_modify_device |= force->bond->datamask_modify;
} else {
datamask_read_device |= force->bond->datamask_read;
datamask_modify_device |= force->bond->datamask_modify;
}
}
if (atomKK->molecular && force->angle) {
if (force->angle->execution_space==Host) {
execute_on_host = true;
datamask_read_host |= force->angle->datamask_read;
datamask_modify_device |= force->angle->datamask_modify;
} else {
datamask_read_device |= force->angle->datamask_read;
datamask_modify_device |= force->angle->datamask_modify;
}
}
if (atomKK->molecular && force->dihedral) {
if (force->dihedral->execution_space==Host) {
execute_on_host = true;
datamask_read_host |= force->dihedral->datamask_read;
datamask_modify_device |= force->dihedral->datamask_modify;
} else {
datamask_read_device |= force->dihedral->datamask_read;
datamask_modify_device |= force->dihedral->datamask_modify;
}
}
if (atomKK->molecular && force->improper) {
if (force->improper->execution_space==Host) {
execute_on_host = true;
datamask_read_host |= force->improper->datamask_read;
datamask_modify_device |= force->improper->datamask_modify;
} else {
datamask_read_device |= force->improper->datamask_read;
datamask_modify_device |= force->improper->datamask_modify;
}
}
if (kspace_compute_flag) {
if (force->kspace->execution_space==Host) {
execute_on_host = true;
datamask_read_host |= force->kspace->datamask_read;
datamask_modify_device |= force->kspace->datamask_modify;
} else {
datamask_read_device |= force->kspace->datamask_read;
datamask_modify_device |= force->kspace->datamask_modify;
}
}
if (pair_compute_flag) {
atomKK->sync(force->pair->execution_space,force->pair->datamask_read);
atomKK->sync(force->pair->execution_space,~(~force->pair->datamask_read|(F_MASK | ENERGY_MASK | VIRIAL_MASK)));
Kokkos::Timer ktimer;
force->pair->compute(eflag,vflag);
atomKK->modified(force->pair->execution_space,force->pair->datamask_modify);
atomKK->modified(force->pair->execution_space,~(~force->pair->datamask_modify|(F_MASK | ENERGY_MASK | VIRIAL_MASK)));
timer->stamp(Timer::PAIR);
}
if (execute_on_host) {
if (pair_compute_flag && force->pair->datamask_modify!=(F_MASK | ENERGY_MASK | VIRIAL_MASK))
Kokkos::fence();
atomKK->sync_overlapping_device(Host,~(~datamask_read_host|(F_MASK | ENERGY_MASK | VIRIAL_MASK)));
if (pair_compute_flag && force->pair->execution_space!=Host) {
Kokkos::deep_copy(LMPHostType(),atomKK->k_f.h_view,0.0);
}
}
if (atomKK->molecular) {
if (force->bond) {
atomKK->sync(force->bond->execution_space,~(~force->bond->datamask_read|(F_MASK | ENERGY_MASK | VIRIAL_MASK)));
force->bond->compute(eflag,vflag);
atomKK->modified(force->bond->execution_space,~(~force->bond->datamask_modify|(F_MASK | ENERGY_MASK | VIRIAL_MASK)));
}
if (force->angle) {
atomKK->sync(force->angle->execution_space,~(~force->angle->datamask_read|(F_MASK | ENERGY_MASK | VIRIAL_MASK)));
force->angle->compute(eflag,vflag);
atomKK->modified(force->angle->execution_space,~(~force->angle->datamask_modify|(F_MASK | ENERGY_MASK | VIRIAL_MASK)));
}
if (force->dihedral) {
atomKK->sync(force->dihedral->execution_space,~(~force->dihedral->datamask_read|(F_MASK | ENERGY_MASK | VIRIAL_MASK)));
force->dihedral->compute(eflag,vflag);
atomKK->modified(force->dihedral->execution_space,~(~force->dihedral->datamask_modify|(F_MASK | ENERGY_MASK | VIRIAL_MASK)));
}
if (force->improper) {
atomKK->sync(force->improper->execution_space,~(~force->improper->datamask_read|(F_MASK | ENERGY_MASK | VIRIAL_MASK)));
force->improper->compute(eflag,vflag);
atomKK->modified(force->improper->execution_space,~(~force->improper->datamask_modify|(F_MASK | ENERGY_MASK | VIRIAL_MASK)));
}
timer->stamp(Timer::BOND);
}
if (kspace_compute_flag) {
atomKK->sync(force->kspace->execution_space,~(~force->kspace->datamask_read|(F_MASK | ENERGY_MASK | VIRIAL_MASK)));
force->kspace->compute(eflag,vflag);
atomKK->modified(force->kspace->execution_space,~(~force->kspace->datamask_modify|(F_MASK | ENERGY_MASK | VIRIAL_MASK)));
timer->stamp(Timer::KSPACE);
}
if (execute_on_host && !std::is_same<LMPHostType,LMPDeviceType>::value) {
if (f_merge_copy.extent(0)<atomKK->k_f.extent(0)) {
f_merge_copy = DAT::t_f_array("DynamicalMatrixKokkos::f_merge_copy",atomKK->k_f.extent(0));
}
f = atomKK->k_f.d_view;
Kokkos::deep_copy(LMPHostType(),f_merge_copy,atomKK->k_f.h_view);
Kokkos::parallel_for(atomKK->k_f.extent(0),
ForceAdder<DAT::t_f_array,DAT::t_f_array>(atomKK->k_f.d_view,f_merge_copy));
atomKK->k_f.clear_sync_state(); // special case
atomKK->k_f.modify<LMPDeviceType>();
}
if (n_pre_reverse) {
modify->pre_reverse(eflag,vflag);
timer->stamp(Timer::MODIFY);
}
if (force->newton) {
comm->reverse_comm();
timer->stamp(Timer::COMM);
}
// force modifications
if (n_post_force) {
modify->post_force(vflag);
timer->stamp(Timer::MODIFY);
}
atomKK->sync(Host,F_MASK);
lmp->kokkos->auto_sync = 1;
++ update->nsteps;
}
/* ----------------------------------------------------------------------
clear force on own & ghost atoms
clear other arrays as needed
------------------------------------------------------------------------- */
void DynamicalMatrixKokkos::force_clear()
{
if (external_force_clear) return;
atomKK->k_f.clear_sync_state(); // ignore host forces/torques since device views
// clear force on all particles
// if either newton flag is set, also include ghosts
// when using threads always clear all forces.
int nall = atomKK->nlocal;
if (force->newton) nall += atomKK->nghost;
Kokkos::parallel_for(nall, Zero<typename ArrayTypes<LMPDeviceType>::t_f_array>(atomKK->k_f.view<LMPDeviceType>()));
atomKK->modified(Device,F_MASK);
}

54
src/KOKKOS/dynamical_matrix_kokkos.h Normal file
View File

@ -0,0 +1,54 @@
/* -*- c++ -*- ----------------------------------------------------------
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.
------------------------------------------------------------------------- */
#ifdef COMMAND_CLASS
// clang-format off
CommandStyle(dynamical_matrix/kk,DynamicalMatrixKokkos);
CommandStyle(dynamical_matrix/kk/device,DynamicalMatrixKokkos);
CommandStyle(dynamical_matrix/kk/host,DynamicalMatrixKokkos);
// clang-format on
#else
#ifndef LMP_DYNAMICAL_MATRIX_KOKKOS_H
#define LMP_DYNAMICAL_MATRIX_KOKKOS_H
#include "dynamical_matrix.h"
#include "kokkos_type.h"
namespace LAMMPS_NS {
class DynamicalMatrixKokkos : public DynamicalMatrix {
public:
DynamicalMatrixKokkos(class LAMMPS *);
~DynamicalMatrixKokkos() override;
void command(int, char **) override;
void setup();
KOKKOS_INLINE_FUNCTION
void operator() (const int& i) const {
f(i,0) += f_merge_copy(i,0);
f(i,1) += f_merge_copy(i,1);
f(i,2) += f_merge_copy(i,2);
}
protected:
void update_force();
void force_clear();
DAT::t_f_array f_merge_copy,f;
};
} // namespace LAMMPS_NS
#endif //LMP_DYNAMICAL_MATRIX_KOKKOS_H
#endif

353
src/KOKKOS/third_order_kokkos.cpp Normal file
View File

@ -0,0 +1,353 @@
// clang-format off
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Charlie Sievers (UC Davis), charliesievers at cox.net
------------------------------------------------------------------------- */
#include "third_order_kokkos.h"
#include "angle.h"
#include "atom_kokkos.h"
#include "atom_masks.h"
#include "bond.h"
#include "comm.h"
#include "compute.h"
#include "dihedral.h"
#include "domain.h"
#include "error.h"
#include "finish.h"
#include "force.h"
#include "group.h"
#include "improper.h"
#include "kokkos.h"
#include "kspace.h"
#include "memory.h"
#include "modify.h"
#include "neighbor.h"
#include "pair.h"
#include "timer.h"
#include "update.h"
#include <cmath>
#include <cstring>
#include <algorithm>
using namespace LAMMPS_NS;
enum{REGULAR,ESKM};
template<class ViewA, class ViewB>
struct ForceAdder {
ViewA a;
ViewB b;
ForceAdder(const ViewA& a_, const ViewB& b_):a(a_),b(b_) {}
KOKKOS_INLINE_FUNCTION
void operator() (const int& i) const {
a(i,0) += b(i,0);
a(i,1) += b(i,1);
a(i,2) += b(i,2);
}
};
/* ---------------------------------------------------------------------- */
template<class View>
struct Zero {
View v;
Zero(const View &v_):v(v_) {}
KOKKOS_INLINE_FUNCTION
void operator()(const int &i) const {
v(i,0) = 0;
v(i,1) = 0;
v(i,2) = 0;
}
};
/* ---------------------------------------------------------------------- */
ThirdOrderKokkos::ThirdOrderKokkos(LAMMPS *lmp) : ThirdOrder(lmp)
{
atomKK = (AtomKokkos *) atom;
}
/* ---------------------------------------------------------------------- */
ThirdOrderKokkos::~ThirdOrderKokkos()
{
}
/* ---------------------------------------------------------------------- */
void ThirdOrderKokkos::command(int narg, char **arg)
{
atomKK->sync(Host, X_MASK|RMASS_MASK|TYPE_MASK);
ThirdOrder::command(narg, arg);
}
/* ----------------------------------------------------------------------
setup without output or one-time post-init setup
flag = 0 = just force calculation
flag = 1 = reneighbor and force calculation
------------------------------------------------------------------------- */
void ThirdOrderKokkos::setup()
{
lmp->kokkos->auto_sync = 1;
// setup domain, communication and neighboring
// acquire ghosts
// build neighbor lists
if (triclinic) domain->x2lamda(atom->nlocal);
domain->pbc();
domain->reset_box();
comm->setup();
if (neighbor->style) neighbor->setup_bins();
comm->exchange();
comm->borders();
if (triclinic) domain->lamda2x(atom->nlocal+atom->nghost);
domain->image_check();
domain->box_too_small_check();
neighbor->build(1);
// compute all forces
eflag=0;
vflag=0;
if (force->kspace) {
force->kspace->setup();
}
update_force();
if (pair_compute_flag) {
atomKK->sync(force->pair->execution_space,force->pair->datamask_read);
force->pair->compute(eflag,vflag);
atomKK->modified(force->pair->execution_space,force->pair->datamask_modify);
}
else if (force->pair) force->pair->compute_dummy(eflag,vflag);
update->setupflag = 0;
lmp->kokkos->auto_sync = 0;
//if all then skip communication groupmap population
if (gcount == atom->natoms)
for (bigint i=0; i<atom->natoms; i++)
groupmap[i] = i;
else
create_groupmap();
}
/* ----------------------------------------------------------------------
evaluate potential energy and forces
may migrate atoms due to reneighboring
return new energy, which should include nextra_global dof
return negative gradient stored in atom->f
return negative gradient for nextra_global dof in fextra
------------------------------------------------------------------------- */
void ThirdOrderKokkos::update_force()
{
int n_pre_force = modify->n_pre_force;
int n_pre_reverse = modify->n_pre_reverse;
int n_post_force = modify->n_post_force;
lmp->kokkos->auto_sync = 0;
f_merge_copy = DAT::t_f_array("ThirdOrderKokkos::f_merge_copy",atomKK->k_f.extent(0));
atomKK->modified(Host,X_MASK);
atomKK->sync(Device,X_MASK);
force_clear();
neighbor->ago = 0;
if ((modify->get_fix_by_id("package_intel")) ? true : false)
neighbor->decide();
if (n_pre_force) {
modify->pre_force(vflag);
timer->stamp(Timer::MODIFY);
}
bool execute_on_host = false;
unsigned int datamask_read_device = 0;
unsigned int datamask_modify_device = 0;
unsigned int datamask_read_host = 0;
if (pair_compute_flag) {
if (force->pair->execution_space==Host) {
execute_on_host = true;
datamask_read_host |= force->pair->datamask_read;
datamask_modify_device |= force->pair->datamask_modify;
} else {
datamask_read_device |= force->pair->datamask_read;
datamask_modify_device |= force->pair->datamask_modify;
}
}
if (atomKK->molecular && force->bond) {
if (force->bond->execution_space==Host) {
execute_on_host = true;
datamask_read_host |= force->bond->datamask_read;
datamask_modify_device |= force->bond->datamask_modify;
} else {
datamask_read_device |= force->bond->datamask_read;
datamask_modify_device |= force->bond->datamask_modify;
}
}
if (atomKK->molecular && force->angle) {
if (force->angle->execution_space==Host) {
execute_on_host = true;
datamask_read_host |= force->angle->datamask_read;
datamask_modify_device |= force->angle->datamask_modify;
} else {
datamask_read_device |= force->angle->datamask_read;
datamask_modify_device |= force->angle->datamask_modify;
}
}
if (atomKK->molecular && force->dihedral) {
if (force->dihedral->execution_space==Host) {
execute_on_host = true;
datamask_read_host |= force->dihedral->datamask_read;
datamask_modify_device |= force->dihedral->datamask_modify;
} else {
datamask_read_device |= force->dihedral->datamask_read;
datamask_modify_device |= force->dihedral->datamask_modify;
}
}
if (atomKK->molecular && force->improper) {
if (force->improper->execution_space==Host) {
execute_on_host = true;
datamask_read_host |= force->improper->datamask_read;
datamask_modify_device |= force->improper->datamask_modify;
} else {
datamask_read_device |= force->improper->datamask_read;
datamask_modify_device |= force->improper->datamask_modify;
}
}
if (kspace_compute_flag) {
if (force->kspace->execution_space==Host) {
execute_on_host = true;
datamask_read_host |= force->kspace->datamask_read;
datamask_modify_device |= force->kspace->datamask_modify;
} else {
datamask_read_device |= force->kspace->datamask_read;
datamask_modify_device |= force->kspace->datamask_modify;
}
}
if (pair_compute_flag) {
atomKK->sync(force->pair->execution_space,force->pair->datamask_read);
atomKK->sync(force->pair->execution_space,~(~force->pair->datamask_read|(F_MASK | ENERGY_MASK | VIRIAL_MASK)));
Kokkos::Timer ktimer;
force->pair->compute(eflag,vflag);
atomKK->modified(force->pair->execution_space,force->pair->datamask_modify);
atomKK->modified(force->pair->execution_space,~(~force->pair->datamask_modify|(F_MASK | ENERGY_MASK | VIRIAL_MASK)));
timer->stamp(Timer::PAIR);
}
if (execute_on_host) {
if (pair_compute_flag && force->pair->datamask_modify!=(F_MASK | ENERGY_MASK | VIRIAL_MASK))
Kokkos::fence();
atomKK->sync_overlapping_device(Host,~(~datamask_read_host|(F_MASK | ENERGY_MASK | VIRIAL_MASK)));
if (pair_compute_flag && force->pair->execution_space!=Host) {
Kokkos::deep_copy(LMPHostType(),atomKK->k_f.h_view,0.0);
}
}
if (atomKK->molecular) {
if (force->bond) {
atomKK->sync(force->bond->execution_space,~(~force->bond->datamask_read|(F_MASK | ENERGY_MASK | VIRIAL_MASK)));
force->bond->compute(eflag,vflag);
atomKK->modified(force->bond->execution_space,~(~force->bond->datamask_modify|(F_MASK | ENERGY_MASK | VIRIAL_MASK)));
}
if (force->angle) {
atomKK->sync(force->angle->execution_space,~(~force->angle->datamask_read|(F_MASK | ENERGY_MASK | VIRIAL_MASK)));
force->angle->compute(eflag,vflag);
atomKK->modified(force->angle->execution_space,~(~force->angle->datamask_modify|(F_MASK | ENERGY_MASK | VIRIAL_MASK)));
}
if (force->dihedral) {
atomKK->sync(force->dihedral->execution_space,~(~force->dihedral->datamask_read|(F_MASK | ENERGY_MASK | VIRIAL_MASK)));
force->dihedral->compute(eflag,vflag);
atomKK->modified(force->dihedral->execution_space,~(~force->dihedral->datamask_modify|(F_MASK | ENERGY_MASK | VIRIAL_MASK)));
}
if (force->improper) {
atomKK->sync(force->improper->execution_space,~(~force->improper->datamask_read|(F_MASK | ENERGY_MASK | VIRIAL_MASK)));
force->improper->compute(eflag,vflag);
atomKK->modified(force->improper->execution_space,~(~force->improper->datamask_modify|(F_MASK | ENERGY_MASK | VIRIAL_MASK)));
}
timer->stamp(Timer::BOND);
}
if (kspace_compute_flag) {
atomKK->sync(force->kspace->execution_space,~(~force->kspace->datamask_read|(F_MASK | ENERGY_MASK | VIRIAL_MASK)));
force->kspace->compute(eflag,vflag);
atomKK->modified(force->kspace->execution_space,~(~force->kspace->datamask_modify|(F_MASK | ENERGY_MASK | VIRIAL_MASK)));
timer->stamp(Timer::KSPACE);
}
if (execute_on_host && !std::is_same<LMPHostType,LMPDeviceType>::value) {
if (f_merge_copy.extent(0)<atomKK->k_f.extent(0)) {
f_merge_copy = DAT::t_f_array("ThirdOrderKokkos::f_merge_copy",atomKK->k_f.extent(0));
}
f = atomKK->k_f.d_view;
Kokkos::deep_copy(LMPHostType(),f_merge_copy,atomKK->k_f.h_view);
Kokkos::parallel_for(atomKK->k_f.extent(0),
ForceAdder<DAT::t_f_array,DAT::t_f_array>(atomKK->k_f.d_view,f_merge_copy));
atomKK->k_f.clear_sync_state(); // special case
atomKK->k_f.modify<LMPDeviceType>();
}
if (n_pre_reverse) {
modify->pre_reverse(eflag,vflag);
timer->stamp(Timer::MODIFY);
}
if (force->newton) {
comm->reverse_comm();
timer->stamp(Timer::COMM);
}
// force modifications
if (n_post_force) {
modify->post_force(vflag);
timer->stamp(Timer::MODIFY);
}
atomKK->sync(Host,F_MASK);
lmp->kokkos->auto_sync = 1;
++ update->nsteps;
}
/* ----------------------------------------------------------------------
clear force on own & ghost atoms
clear other arrays as needed
------------------------------------------------------------------------- */
void ThirdOrderKokkos::force_clear()
{
if (external_force_clear) return;
atomKK->k_f.clear_sync_state(); // ignore host forces/torques since device views
// clear force on all particles
// if either newton flag is set, also include ghosts
// when using threads always clear all forces.
int nall = atomKK->nlocal;
if (force->newton) nall += atomKK->nghost;
Kokkos::parallel_for(nall, Zero<typename ArrayTypes<LMPDeviceType>::t_f_array>(atomKK->k_f.view<LMPDeviceType>()));
atomKK->modified(Device,F_MASK);
}

54
src/KOKKOS/third_order_kokkos.h Normal file
View File

@ -0,0 +1,54 @@
/* -*- c++ -*- ----------------------------------------------------------
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.
------------------------------------------------------------------------- */
#ifdef COMMAND_CLASS
// clang-format off
CommandStyle(third_order/kk,ThirdOrderKokkos);
CommandStyle(third_order/kk/device,ThirdOrderKokkos);
CommandStyle(third_order/kk/host,ThirdOrderKokkos);
// clang-format on
#else
#ifndef LMP_THIRD_ORDER_KOKKOS_H
#define LMP_THIRD_ORDER_KOKKOS_H
#include "third_order.h"
#include "kokkos_type.h"
namespace LAMMPS_NS {
class ThirdOrderKokkos : public ThirdOrder {
public:
ThirdOrderKokkos(class LAMMPS *);
~ThirdOrderKokkos() override;
void command(int, char **) override;
void setup();
KOKKOS_INLINE_FUNCTION
void operator() (const int& i) const {
f(i,0) += f_merge_copy(i,0);
f(i,1) += f_merge_copy(i,1);
f(i,2) += f_merge_copy(i,2);
}
protected:
void update_force();
void force_clear();
DAT::t_f_array f_merge_copy,f;
};
} // namespace LAMMPS_NS
#endif //LMP_THIRD_ORDER_KOKKOS_H
#endif

740
src/PHONON/dynamical_matrix.cpp
View File

@ -10,11 +10,11 @@
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
----------------------------------------------------------------------- */
------------------------------------------------------------------------- */
//
// Created by charlie sievers on 6/21/18.
//
/* ----------------------------------------------------------------------
Contributing author: Charlie Sievers (UC Davis), charliesievers at cox.net
------------------------------------------------------------------------- */
#include "dynamical_matrix.h"
@ -22,6 +22,7 @@
#include "atom.h"
#include "bond.h"
#include "comm.h"
#include "compute.h"
#include "dihedral.h"
#include "domain.h"
#include "error.h"
@ -48,14 +49,14 @@ enum{REGULAR,ESKM};
DynamicalMatrix::DynamicalMatrix(LAMMPS *lmp) : Command(lmp), fp(nullptr)
{
external_force_clear = 1;
external_force_clear = 0;
}
/* ---------------------------------------------------------------------- */
DynamicalMatrix::~DynamicalMatrix()
{
if (fp && me == 0) {
if (fp && comm->me == 0) {
if (compressed) platform::pclose(fp);
else fclose(fp);
memory->destroy(groupmap);
@ -71,111 +72,120 @@ DynamicalMatrix::~DynamicalMatrix()
void DynamicalMatrix::setup()
{
// setup domain, communication and neighboring
// acquire ghosts
// build neighbor lists
if (triclinic) domain->x2lamda(atom->nlocal);
domain->pbc();
domain->reset_box();
comm->setup();
if (neighbor->style) neighbor->setup_bins();
comm->exchange();
comm->borders();
if (triclinic) domain->lamda2x(atom->nlocal+atom->nghost);
domain->image_check();
domain->box_too_small_check();
neighbor->build(1);
// setup domain, communication and neighboring
// acquire ghosts
// build neighbor lists
if (triclinic) domain->x2lamda(atom->nlocal);
domain->pbc();
domain->reset_box();
comm->setup();
if (neighbor->style) neighbor->setup_bins();
comm->exchange();
comm->borders();
if (triclinic) domain->lamda2x(atom->nlocal+atom->nghost);
domain->image_check();
domain->box_too_small_check();
neighbor->build(1);
// compute all forces
external_force_clear = 0;
eflag=0;
vflag=0;
update_force();
// compute all forces
eflag=0;
vflag=0;
if (force->kspace) {
force->kspace->setup();
}
update_force();
//if all then skip communication groupmap population
if (gcount == atom->natoms)
for (bigint i=0; i<atom->natoms; i++)
groupmap[i] = i;
else
create_groupmap();
update->setupflag = 0;
//if all then skip communication groupmap population
if (gcount == atom->natoms)
for (bigint i=0; i<atom->natoms; i++)
groupmap[i] = i;
else
create_groupmap();
}
/* ---------------------------------------------------------------------- */
void DynamicalMatrix::command(int narg, char **arg)
{
MPI_Comm_rank(world,&me);
MPI_Comm_rank(world,&me);
if (domain->box_exist == 0)
error->all(FLERR,"Dynamical_matrix command before simulation box is defined");
if (narg < 2) error->all(FLERR,"Illegal dynamical_matrix command");
if (domain->box_exist == 0)
error->all(FLERR,"Dynamical_matrix command before simulation box is defined");
if (narg < 2) error->all(FLERR,"Illegal dynamical_matrix command");
lmp->init();
lmp->init();
// orthogonal vs triclinic simulation box
// orthogonal vs triclinic simulation box
triclinic = domain->triclinic;
triclinic = domain->triclinic;
if (force->pair && force->pair->compute_flag) pair_compute_flag = 1;
else pair_compute_flag = 0;
if (force->kspace && force->kspace->compute_flag) kspace_compute_flag = 1;
else kspace_compute_flag = 0;
if (force->pair && force->pair->compute_flag) pair_compute_flag = 1;
else pair_compute_flag = 0;
if (force->kspace && force->kspace->compute_flag) kspace_compute_flag = 1;
else kspace_compute_flag = 0;
// group and style
// group and style
igroup = group->find(arg[0]);
if (igroup == -1) error->all(FLERR,"Could not find dynamical matrix group ID");
groupbit = group->bitmask[igroup];
gcount = group->count(igroup);
dynlen = (gcount)*3;
memory->create(groupmap,atom->natoms,"total_group_map:totalgm");
update->setupflag = 1;
igroup = group->find(arg[0]);
if (igroup == -1) error->all(FLERR,"Could not find dynamical matrix group ID");
groupbit = group->bitmask[igroup];
gcount = group->count(igroup);
dynlen = (gcount)*3;
memory->create(groupmap,atom->natoms,"total_group_map:totalgm");
update->setupflag = 1;
int style = -1;
if (strcmp(arg[1],"regular") == 0) style = REGULAR;
else if (strcmp(arg[1],"eskm") == 0) style = ESKM;
else error->all(FLERR,"Illegal Dynamical Matrix command");
del = utils::numeric(FLERR, arg[2],false,lmp);
int style = -1;
if (strcmp(arg[1],"regular") == 0) style = REGULAR;
else if (strcmp(arg[1],"eskm") == 0) style = ESKM;
else error->all(FLERR,"Illegal Dynamical_matrix command");
del = utils::numeric(FLERR, arg[2],false,lmp);
// set option defaults
// set option defaults
binaryflag = 0;
scaleflag = 0;
compressed = 0;
file_flag = 0;
file_opened = 0;
conversion = 1;
binaryflag = 0;
scaleflag = 0;
compressed = 0;
file_flag = 0;
file_opened = 0;
folded = 0;
conversion = 1;
// read options from end of input line
if (style == REGULAR) options(narg-3,&arg[3]); //COME BACK
else if (style == ESKM) options(narg-3,&arg[3]); //COME BACK
else if (comm->me == 0 && screen) fprintf(screen,"Illegal Dynamical Matrix command\n");
// read options from end of input line
if (style == REGULAR) options(narg-3,&arg[3]);
else if (style == ESKM) options(narg-3,&arg[3]);
else if (me == 0 && screen) fprintf(screen,"Illegal Dynamical Matrix command\n");
if (atom->map_style == Atom::MAP_NONE)
error->all(FLERR,"Dynamical_matrix command requires an atom map, see atom_modify");
if (!folded) dynlenb = dynlen;
else dynlenb = (atom->natoms)*3;
// move atoms by 3-vector or specified variable(s)
if (atom->map_style == Atom::MAP_NONE)
error->all(FLERR,"Dynamical_matrix command requires an atom map");
if (style == REGULAR) {
setup();
timer->init();
timer->barrier_start();
calculateMatrix();
timer->barrier_stop();
}
// move atoms by 3-vector or specified variable(s)
if (style == ESKM) {
setup();
convert_units(update->unit_style);
conversion = conv_energy/conv_distance/conv_mass;
timer->init();
timer->barrier_start();
calculateMatrix();
timer->barrier_stop();
}
if (style == REGULAR) {
setup();
timer->init();
timer->barrier_start();
calculateMatrix();
timer->barrier_stop();
}
Finish finish(lmp);
finish.end(1);
if (style == ESKM) {
setup();
convert_units(update->unit_style);
conversion = conv_energy/conv_distance/conv_mass;
timer->init();
timer->barrier_start();
calculateMatrix();
timer->barrier_stop();
}
Finish finish(lmp);
finish.end(1);
}
/* ----------------------------------------------------------------------
@ -192,9 +202,9 @@ void DynamicalMatrix::options(int narg, char **arg)
if (strcmp(arg[iarg],"binary") == 0) {
if (iarg + 2 > narg) error->all(FLERR, "Illegal dynamical_matrix command");
if (strcmp(arg[iarg+1],"gzip") == 0) {
compressed = 1;
compressed = 1;
} else {
binaryflag = utils::logical(FLERR,arg[iarg+1],false,lmp);
binaryflag = utils::logical(FLERR,arg[iarg+1],false,lmp);
}
iarg += 2;
} else if (strcmp(arg[iarg],"file") == 0) {
@ -202,6 +212,14 @@ void DynamicalMatrix::options(int narg, char **arg)
filename = arg[iarg + 1];
file_flag = 1;
iarg += 2;
} else if (strcmp(arg[iarg],"fold") == 0) {
if (iarg+2 > narg) error->all(FLERR, "Illegal dynamical_matrix command");
if (strcmp(arg[iarg+1],"yes") == 0) {
folded = 1;
} else if (strcmp(arg[iarg+1],"no") == 0) {
folded = 0;
} else error->all(FLERR,"Illegal input for dynamical_matrix fold option");
iarg += 2;
} else error->all(FLERR,"Illegal dynamical_matrix command");
}
if (file_flag == 1) {
@ -217,23 +235,23 @@ void DynamicalMatrix::options(int narg, char **arg)
void DynamicalMatrix::openfile(const char *filename)
{
// if file already opened, return
if (file_opened) return;
fp = nullptr;
// if file already opened, return
if (file_opened) return;
fp = nullptr;
if (me == 0) {
if (compressed) {
fp = platform::compressed_write(std::string(filename)+".gz");
if (!fp) error->one(FLERR,"Cannot open compressed file");
} else if (binaryflag) {
fp = fopen(filename,"wb");
} else {
fp = fopen(filename,"w");
}
if (!fp) error->one(FLERR,"Cannot open dynmat file: {}", utils::getsyserror());
if (me == 0) {
if (compressed) {
fp = platform::compressed_write(std::string(filename)+".gz");
if (!fp) error->one(FLERR,"Cannot open compressed file");
} else if (binaryflag) {
fp = fopen(filename,"wb");
} else {
fp = fopen(filename,"w");
}
if (!fp) error->one(FLERR,"Cannot open dynmat file: {}", utils::getsyserror());
}
file_opened = 1;
file_opened = 1;
}
/* ----------------------------------------------------------------------
@ -242,98 +260,112 @@ void DynamicalMatrix::openfile(const char *filename)
void DynamicalMatrix::calculateMatrix()
{
int local_idx; // local index
int local_jdx; // second local index
int nlocal = atom->nlocal;
bigint natoms = atom->natoms;
int *type = atom->type;
bigint *gm = groupmap;
double imass; // dynamical matrix element
double *m = atom->mass;
double **f = atom->f;
int local_idx; // local index
int local_jdx; // second local index
int nlocal = atom->nlocal;
bigint natoms = atom->natoms;
int *type = atom->type;
bigint *gm = groupmap;
double imass; // dynamical matrix element
double *m = atom->mass;
double **f = atom->f;
double **dynmat = new double*[3];
for (int i=0; i<3; i++)
dynmat[i] = new double[dynlen];
double **dynmat = new double*[3];
for (int i=0; i<3; i++)
dynmat[i] = new double[dynlenb];
double **fdynmat = new double*[3];
for (int i=0; i<3; i++)
fdynmat[i] = new double[dynlen];
double **fdynmat = new double*[3];
for (int i=0; i<3; i++)
fdynmat[i] = new double[dynlenb];
//initialize dynmat to all zeros
//initialize dynmat to all zeros
dynmat_clear(dynmat);
if (me == 0 && screen) {
fprintf(screen,"Calculating Dynamical Matrix ...\n");
fprintf(screen," Total # of atoms = " BIGINT_FORMAT "\n", natoms);
fprintf(screen," Atoms in group = " BIGINT_FORMAT "\n", gcount);
fprintf(screen," Total dynamical matrix elements = " BIGINT_FORMAT "\n", (dynlen*dynlen) );
}
// emit dynlen rows of dimalpha*dynlen*dimbeta elements
update->nsteps = 0;
int prog = 0;
for (bigint i=1; i<=natoms; i++) {
local_idx = atom->map(i);
if (gm[i-1] < 0)
continue;
for (int alpha=0; alpha<3; alpha++) {
displace_atom(local_idx, alpha, 1);
update_force();
for (bigint j=1; j<=natoms; j++) {
local_jdx = atom->map(j);
if (local_idx >= 0 && local_jdx >= 0 && local_jdx < nlocal
&& (gm[j-1] >= 0 || folded)){
if (folded) {
for (int beta=0; beta<3; beta++){
dynmat[alpha][(j-1)*3+beta] -= f[local_jdx][beta];
}
} else {
for (int beta=0; beta<3; beta++){
dynmat[alpha][gm[j-1]*3+beta] -= f[local_jdx][beta];
}
}
}
}
displace_atom(local_idx,alpha,-2);
update_force();
for (bigint j=1; j<=natoms; j++) {
local_jdx = atom->map(j);
if (local_idx >= 0 && local_jdx >= 0 && local_jdx < nlocal
&& (gm[j-1] >= 0 || folded)){
if (atom->rmass_flag == 1)
imass = sqrt(m[local_idx] * m[local_jdx]);
else
imass = sqrt(m[type[local_idx]] * m[type[local_jdx]]);
if (folded){
for (int beta=0; beta<3; beta++){
dynmat[alpha][(j-1)*3+beta] -= -f[local_jdx][beta];
dynmat[alpha][(j-1)*3+beta] /= (2 * del * imass);
dynmat[alpha][(j-1)*3+beta] *= conversion;
}
} else {
for (int beta=0; beta<3; beta++){
dynmat[alpha][gm[j-1]*3+beta] -= -f[local_jdx][beta];
dynmat[alpha][gm[j-1]*3+beta] /= (2 * del * imass);
dynmat[alpha][gm[j-1]*3+beta] *= conversion;
}
}
}
}
displace_atom(local_idx,alpha,1);
}
for (int k=0; k<3; k++)
MPI_Reduce(dynmat[k],fdynmat[k],dynlenb,MPI_DOUBLE,MPI_SUM,0,world);
if (me == 0)
writeMatrix(fdynmat);
dynmat_clear(dynmat);
if (comm->me == 0 && screen) {
fprintf(screen,"Calculating Dynamical Matrix ...\n");
fprintf(screen," Total # of atoms = " BIGINT_FORMAT "\n", natoms);
fprintf(screen," Atoms in group = " BIGINT_FORMAT "\n", gcount);
fprintf(screen," Total dynamical matrix elements = " BIGINT_FORMAT "\n", (dynlen*dynlen) );
if (me == 0 && screen) {
int p = 10 * gm[i-1] / gcount;
if (p > prog) {
prog = p;
fprintf(screen," %d%%",p*10);
fflush(screen);
}
}
}
if (me == 0 && screen) fprintf(screen,"\n");
// emit dynlen rows of dimalpha*dynlen*dimbeta elements
for (int i=0; i < 3; i++)
delete [] dynmat[i];
delete [] dynmat;
update->nsteps = 0;
int prog = 0;
for (bigint i=1; i<=natoms; i++) {
local_idx = atom->map(i);
if (gm[i-1] < 0)
continue;
for (int alpha=0; alpha<3; alpha++) {
displace_atom(local_idx, alpha, 1);
update_force();
for (bigint j=1; j<=natoms; j++) {
local_jdx = atom->map(j);
if (local_idx >= 0 && local_jdx >= 0 && local_jdx < nlocal
&& gm[j-1] >= 0) {
for (int beta=0; beta<3; beta++) {
dynmat[alpha][gm[j-1]*3+beta] -= f[local_jdx][beta];
}
}
}
displace_atom(local_idx,alpha,-2);
update_force();
for (bigint j=1; j<=natoms; j++) {
local_jdx = atom->map(j);
if (local_idx >= 0 && local_jdx >= 0 && local_jdx < nlocal
&& gm[j-1] >= 0) {
for (int beta=0; beta<3; beta++) {
if (atom->rmass_flag == 1)
imass = sqrt(m[local_idx] * m[local_jdx]);
else
imass = sqrt(m[type[local_idx]] * m[type[local_jdx]]);
dynmat[alpha][gm[j-1]*3+beta] -= -f[local_jdx][beta];
dynmat[alpha][gm[j-1]*3+beta] /= (2 * del * imass);
dynmat[alpha][gm[j-1]*3+beta] *= conversion;
}
}
}
displace_atom(local_idx,alpha,1);
}
for (int k=0; k<3; k++)
MPI_Reduce(dynmat[k],fdynmat[k],dynlen,MPI_DOUBLE,MPI_SUM,0,world);
if (me == 0)
writeMatrix(fdynmat);
dynmat_clear(dynmat);
if (comm->me == 0 && screen) {
int p = 10 * gm[i-1] / gcount;
if (p > prog) {
prog = p;
fprintf(screen," %d%%",p*10);
fflush(screen);
}
}
}
if (comm->me == 0 && screen) fprintf(screen,"\n");
for (int i=0; i < 3; i++)
delete [] fdynmat[i];
delete [] fdynmat;
for (int i=0; i < 3; i++)
delete [] dynmat[i];
delete [] dynmat;
for (int i=0; i < 3; i++)
delete [] fdynmat[i];
delete [] fdynmat;
if (screen && me ==0 ) fprintf(screen,"Finished Calculating Dynamical Matrix\n");
if (screen && me ==0 ) fprintf(screen,"Finished Calculating Dynamical Matrix\n");
}
/* ----------------------------------------------------------------------
@ -342,25 +374,25 @@ void DynamicalMatrix::calculateMatrix()
void DynamicalMatrix::writeMatrix(double **dynmat)
{
if (me != 0 || !fp)
return;
if (me != 0 || !fp)
return;
clearerr(fp);
if (binaryflag) {
for (int i=0; i<3; i++)
fwrite(dynmat[i], sizeof(double), dynlen, fp);
if (ferror(fp))
error->one(FLERR, "Error writing to binary file");
} else {
for (int i = 0; i < 3; i++) {
for (bigint j = 0; j < dynlen; j++) {
if ((j+1)%3==0) fprintf(fp, "%4.8f\n", dynmat[i][j]);
else fprintf(fp, "%4.8f ",dynmat[i][j]);
}
}
if (ferror(fp))
error->one(FLERR,"Error writing to file");
clearerr(fp);
if (binaryflag) {
for (int i=0; i<3; i++)
fwrite(dynmat[i], sizeof(double), dynlenb, fp);
if (ferror(fp))
error->one(FLERR, "Error writing to binary file");
} else {
for (int i = 0; i < 3; i++) {
for (bigint j = 0; j < dynlenb; j++) {
if ((j+1)%3==0) fprintf(fp, "%4.8f\n", dynmat[i][j]);
else fprintf(fp, "%4.8f ",dynmat[i][j]);
}
}
if (ferror(fp))
error->one(FLERR,"Error writing to file");
}
}
/* ----------------------------------------------------------------------
@ -369,18 +401,18 @@ void DynamicalMatrix::writeMatrix(double **dynmat)
void DynamicalMatrix::displace_atom(int local_idx, int direction, int magnitude)
{
if (local_idx < 0) return;
if (local_idx < 0) return;
double **x = atom->x;
int *sametag = atom->sametag;
int j = local_idx;
double **x = atom->x;
int *sametag = atom->sametag;
int j = local_idx;
x[local_idx][direction] += del*magnitude;
x[local_idx][direction] += del*magnitude;
while (sametag[j] >= 0) {
j = sametag[j];
x[j][direction] += del*magnitude;
}
while (sametag[j] >= 0) {
j = sametag[j];
x[j][direction] += del*magnitude;
}
}
@ -394,35 +426,50 @@ void DynamicalMatrix::displace_atom(int local_idx, int direction, int magnitude)
void DynamicalMatrix::update_force()
{
force_clear();
int n_post_force = modify->n_post_force;
neighbor->ago = 0;
if ((modify->get_fix_by_id("package_intel")) ? true : false)
neighbor->decide();
force_clear();
int n_pre_force = modify->n_pre_force;
int n_pre_reverse = modify->n_pre_reverse;
int n_post_force = modify->n_post_force;
if (pair_compute_flag) {
force->pair->compute(eflag,vflag);
timer->stamp(Timer::PAIR);
}
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);
timer->stamp(Timer::BOND);
}
if (kspace_compute_flag) {
force->kspace->compute(eflag,vflag);
timer->stamp(Timer::KSPACE);
}
if (force->newton) {
comm->reverse_comm();
timer->stamp(Timer::COMM);
}
// force modifications
if (n_post_force) modify->post_force(vflag);
if (n_pre_force) {
modify->pre_force(vflag);
timer->stamp(Timer::MODIFY);
}
++ update->nsteps;
if (pair_compute_flag) {
force->pair->compute(eflag,vflag);
timer->stamp(Timer::PAIR);
}
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);
timer->stamp(Timer::BOND);
}
if (kspace_compute_flag) {
force->kspace->compute(eflag,vflag);
timer->stamp(Timer::KSPACE);
}
if (n_pre_reverse) {
modify->pre_reverse(eflag,vflag);
timer->stamp(Timer::MODIFY);
}
if (force->newton) {
comm->reverse_comm();
timer->stamp(Timer::COMM);
}
// force modifications
if (n_post_force) {
modify->post_force(vflag);
timer->stamp(Timer::MODIFY);
}
++ update->nsteps;
}
/* ----------------------------------------------------------------------
@ -432,17 +479,17 @@ void DynamicalMatrix::update_force()
void DynamicalMatrix::force_clear()
{
if (external_force_clear) return;
if (external_force_clear) return;
// clear global force array
// if either newton flag is set, also include ghosts
// clear global force array
// if either newton flag is set, also include ghosts
size_t nbytes = sizeof(double) * atom->nlocal;
if (force->newton) nbytes += sizeof(double) * atom->nghost;
size_t nbytes = sizeof(double) * atom->nlocal;
if (force->newton) nbytes += sizeof(double) * atom->nghost;
if (nbytes) {
memset(&atom->f[0][0],0,3*nbytes);
}
if (nbytes) {
memset(&atom->f[0][0],0,3*nbytes);
}
}
/* ----------------------------------------------------------------------
@ -451,67 +498,66 @@ void DynamicalMatrix::force_clear()
void DynamicalMatrix::dynmat_clear(double **dynmat)
{
size_t nbytes = sizeof(double) * dynlenb;
size_t nbytes = sizeof(double) * dynlen;
if (nbytes) {
for (int i=0; i<3; i++)
memset(&dynmat[i][0],0,nbytes);
}
if (nbytes) {
for (int i=0; i<3; i++)
memset(&dynmat[i][0],0,nbytes);
}
}
/* ---------------------------------------------------------------------- */
void DynamicalMatrix::convert_units(const char *style)
{
// physical constants from:
// https://physics.nist.gov/cuu/Constants/Table/allascii.txt
// using thermochemical calorie = 4.184 J
// physical constants from:
// https://physics.nist.gov/cuu/Constants/Table/allascii.txt
// using thermochemical calorie = 4.184 J
if (strcmp(style,"lj") == 0) {
error->all(FLERR,"Conversion Not Set");
//conversion = 1; // lj -> 10 J/mol
if (strcmp(style,"lj") == 0) {
error->all(FLERR,"Conversion Not Set");
//conversion = 1; // lj -> 10 J/mol
} else if (strcmp(style,"real") == 0) {
conv_energy = 418.4; // kcal/mol -> 10 J/mol
conv_mass = 1; // g/mol -> g/mol
conv_distance = 1; // angstrom -> angstrom
} else if (strcmp(style,"real") == 0) {
conv_energy = 418.4; // kcal/mol -> 10 J/mol
conv_mass = 1; // g/mol -> g/mol
conv_distance = 1; // angstrom -> angstrom
} else if (strcmp(style,"metal") == 0) {
conv_energy = 9648.5; // eV -> 10 J/mol
conv_mass = 1; // g/mol -> g/mol
conv_distance = 1; // angstrom -> angstrom
} else if (strcmp(style,"metal") == 0) {
conv_energy = 9648.5; // eV -> 10 J/mol
conv_mass = 1; // g/mol -> g/mol
conv_distance = 1; // angstrom -> angstrom
} else if (strcmp(style,"si") == 0) {
if (comm->me) error->warning(FLERR,"Conversion Warning: Multiplication by Large Float");
conv_energy = 6.022E22; // J -> 10 J/mol
conv_mass = 6.022E26; // kg -> g/mol
conv_distance = 1E-10; // meter -> angstrom
} else if (strcmp(style,"si") == 0) {
if (me) error->warning(FLERR,"Conversion Warning: Multiplication by Large Float");
conv_energy = 6.022E22; // J -> 10 J/mol
conv_mass = 6.022E26; // kg -> g/mol
conv_distance = 1E-10; // meter -> angstrom
} else if (strcmp(style,"cgs") == 0) {
if (comm->me) error->warning(FLERR,"Conversion Warning: Multiplication by Large Float");
conv_energy = 6.022E12; // Erg -> 10 J/mol
conv_mass = 6.022E23; // g -> g/mol
conv_distance = 1E-7; // centimeter -> angstrom
} else if (strcmp(style,"cgs") == 0) {
if (me) error->warning(FLERR,"Conversion Warning: Multiplication by Large Float");
conv_energy = 6.022E12; // Erg -> 10 J/mol
conv_mass = 6.022E23; // g -> g/mol
conv_distance = 1E-7; // centimeter -> angstrom
} else if (strcmp(style,"electron") == 0) {
conv_energy = 262550; // Hartree -> 10 J/mol
conv_mass = 1; // amu -> g/mol
conv_distance = 0.529177249; // bohr -> angstrom
} else if (strcmp(style,"electron") == 0) {
conv_energy = 262550; // Hartree -> 10 J/mol
conv_mass = 1; // amu -> g/mol
conv_distance = 0.529177249; // bohr -> angstrom
} else if (strcmp(style,"micro") == 0) {
if (comm->me) error->warning(FLERR,"Conversion Warning: Untested Conversion");
conv_energy = 6.022E10; // picogram-micrometer^2/microsecond^2 -> 10 J/mol
conv_mass = 6.022E11; // pg -> g/mol
conv_distance = 1E-4; // micrometer -> angstrom
} else if (strcmp(style,"micro") == 0) {
if (me) error->warning(FLERR,"Conversion Warning: Untested Conversion");
conv_energy = 6.022E10; // picogram-micrometer^2/microsecond^2 -> 10 J/mol
conv_mass = 6.022E11; // pg -> g/mol
conv_distance = 1E-4; // micrometer -> angstrom
} else if (strcmp(style,"nano") == 0) {
if (comm->me) error->warning(FLERR,"Conversion Warning: Untested Conversion");
conv_energy = 6.022E4; // attogram-nanometer^2/nanosecond^2 -> 10 J/mol
conv_mass = 6.022E5; // ag -> g/mol
conv_distance = 0.1; // angstrom -> angstrom
} else if (strcmp(style,"nano") == 0) {
if (me) error->warning(FLERR,"Conversion Warning: Untested Conversion");
conv_energy = 6.022E4; // attogram-nanometer^2/nanosecond^2 -> 10 J/mol
conv_mass = 6.022E5; // ag -> g/mol
conv_distance = 0.1; // angstrom -> angstrom
} else error->all(FLERR,"Units Type Conversion Not Found");
} else error->all(FLERR,"Units Type Conversion Not Found");
}
@ -519,66 +565,66 @@ void DynamicalMatrix::convert_units(const char *style)
void DynamicalMatrix::create_groupmap()
{
//Create a group map which maps atom order onto group
// groupmap[global atom index-1] = output column/row
//Create a group map which maps atom order onto group
// groupmap[global atom index-1] = output column/row
int local_idx; // local index
int gid = 0; //group index
int nlocal = atom->nlocal;
int *mask = atom->mask;
bigint natoms = atom->natoms;
int *recv = new int[comm->nprocs];
int *displs = new int[comm->nprocs];
bigint *temp_groupmap = new bigint[natoms];
int local_idx; // local index
int gid = 0; //group index
int nlocal = atom->nlocal;
int *mask = atom->mask;
bigint natoms = atom->natoms;
int *recv = new int[comm->nprocs];
int *displs = new int[comm->nprocs];
bigint *temp_groupmap = new bigint[natoms];
//find number of local atoms in the group (final_gid)
for (bigint i=1; i<=natoms; i++) {
local_idx = atom->map(i);
if ((local_idx >= 0) && (local_idx < nlocal) && mask[local_idx] & groupbit)
gid += 1; // gid at the end of loop is final_Gid
//find number of local atoms in the group (final_gid)
for (bigint i=1; i<=natoms; i++) {
local_idx = atom->map(i);
if ((local_idx >= 0) && (local_idx < nlocal) && mask[local_idx] & groupbit)
gid += 1; // gid at the end of loop is final_Gid
}
//create an array of length final_gid
bigint *sub_groupmap = new bigint[gid];
gid = 0;
//create a map between global atom id and group atom id for each proc
for (bigint i=1; i<=natoms; i++) {
local_idx = atom->map(i);
if ((local_idx >= 0) && (local_idx < nlocal) && mask[local_idx] & groupbit) {
sub_groupmap[gid] = i;
gid += 1;
}
//create an array of length final_gid
bigint *sub_groupmap = new bigint[gid];
}
gid = 0;
//create a map between global atom id and group atom id for each proc
for (bigint i=1; i<=natoms; i++) {
local_idx = atom->map(i);
if ((local_idx >= 0) && (local_idx < nlocal) && mask[local_idx] & groupbit) {
sub_groupmap[gid] = i;
gid += 1;
}
}
//populate arrays for Allgatherv
for (int i=0; i < comm->nprocs; i++) {
recv[i] = 0;
}
recv[me] = gid;
MPI_Allreduce(recv,displs,comm->nprocs,MPI_INT,MPI_SUM,world);
for (int i=0; i<comm->nprocs; i++) {
recv[i]=displs[i];
if (i>0) displs[i] = displs[i-1]+recv[i-1];
else displs[i] = 0;
}
//populate arrays for Allgatherv
for (int i=0; i<comm->nprocs; i++) {
recv[i] = 0;
}
recv[comm->me] = gid;
MPI_Allreduce(recv,displs,comm->nprocs,MPI_INT,MPI_SUM,world);
for (int i=0; i<comm->nprocs; i++) {
recv[i]=displs[i];
if (i>0) displs[i] = displs[i-1]+recv[i-1];
else displs[i] = 0;
}
//combine subgroup maps into total temporary groupmap
MPI_Allgatherv(sub_groupmap,gid,MPI_LMP_BIGINT,temp_groupmap,recv,displs,MPI_LMP_BIGINT,world);
std::sort(temp_groupmap,temp_groupmap+gcount);
//combine subgroup maps into total temporary groupmap
MPI_Allgatherv(sub_groupmap,gid,MPI_LMP_BIGINT,temp_groupmap,recv,displs,MPI_LMP_BIGINT,world);
std::sort(temp_groupmap,temp_groupmap+gcount);
//populate member groupmap based on temp groupmap
bigint j = 0;
for (bigint i=1; i<=natoms; i++) {
// flag groupmap contents that are in temp_groupmap
if (j < gcount && i == temp_groupmap[j])
groupmap[i-1] = j++;
else
groupmap[i-1] = -1;
}
//populate member groupmap based on temp groupmap
bigint j = 0;
for (bigint i=1; i<=natoms; i++) {
// flag groupmap contents that are in temp_groupmap
if (j < gcount && i == temp_groupmap[j])
groupmap[i-1] = j++;
else
groupmap[i-1] = -1;
}
//free that memory!
delete[] recv;
delete[] displs;
delete[] sub_groupmap;
delete[] temp_groupmap;
//free that memory!
delete[] recv;
delete[] displs;
delete[] sub_groupmap;
delete[] temp_groupmap;
}

7
src/PHONON/dynamical_matrix.h
View File

@ -34,11 +34,10 @@ class DynamicalMatrix : public Command {
int nvec; // local atomic dof = length of xvec
void update_force();
void force_clear();
virtual void update_force();
virtual void force_clear();
virtual void openfile(const char *filename);
private:
void options(int, char **);
void calculateMatrix();
void dynmat_clear(double **dynmat);
@ -55,6 +54,7 @@ class DynamicalMatrix : public Command {
int igroup, groupbit;
bigint gcount; // number of atoms in group
bigint dynlen; // rank of dynamical matrix
bigint dynlenb; // new dynlen if folded
int scaleflag;
int me;
bigint *groupmap;
@ -63,6 +63,7 @@ class DynamicalMatrix : public Command {
int binaryflag; // 1 if dump file is written binary, 0 no
int file_opened; // 1 if openfile method has been called, 0 no
int file_flag; // 1 custom file name, 0 dynmat.dat
int folded; // 1 folded, 0 nonfolded
FILE *fp;
};

382
src/PHONON/third_order.cpp
View File

@ -10,11 +10,11 @@
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
----------------------------------------------------------------------- */
------------------------------------------------------------------------- */
//
// Created by charlie sievers on 7/5/18.
//
/* ----------------------------------------------------------------------
Contributing author: Charlie Sievers (UC Davis), charliesievers at cox.net
------------------------------------------------------------------------- */
#include "third_order.h"
@ -32,6 +32,9 @@
#include "kspace.h"
#include "math_special.h"
#include "memory.h"
#include "modify.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "neighbor.h"
#include "pair.h"
#include "timer.h"
@ -42,7 +45,7 @@
using namespace LAMMPS_NS;
using namespace MathSpecial;
enum{REGULAR,BALLISTICO};
enum{REGULAR,ESKM};
/* ---------------------------------------------------------------------- */
@ -86,12 +89,22 @@ void ThirdOrder::setup()
domain->box_too_small_check();
neighbor->build(1);
// build neighbor list this command needs based on earlier request
neighbor->build_one(list);
// compute all forces
external_force_clear = 0;
eflag=0;
vflag=0;
if (force->kspace) {
force->kspace->setup();
}
update_force();
modify->setup(vflag);
update->setupflag = 0;
if (gcount == atom->natoms)
for (bigint i=0; i<atom->natoms; i++)
groupmap[i] = i;
@ -109,7 +122,18 @@ void ThirdOrder::command(int narg, char **arg)
error->all(FLERR,"third_order command before simulation box is defined");
if (narg < 2) error->all(FLERR,"Illegal third_order command");
// request a full neighbor list for use by this command
int irequest = neighbor->request(this);
neighbor->requests[irequest]->pair = 0;
neighbor->requests[irequest]->command = 1;
neighbor->requests[irequest]->half = 0;
neighbor->requests[irequest]->full = 1;
neighbor->requests[irequest]->occasional = 1;
neighbor->requests[irequest]->command_style = "third_order";
lmp->init();
list = neighbor->lists[irequest];
// orthogonal vs triclinic simulation box
@ -123,7 +147,7 @@ void ThirdOrder::command(int narg, char **arg)
// group and style
igroup = group->find(arg[0]);
if (igroup == -1) error->all(FLERR,"Could not find dynamical matrix group ID");
if (igroup == -1) error->all(FLERR,"Could not find third_order group ID");
groupbit = group->bitmask[igroup];
gcount = group->count(igroup);
dynlen = (gcount)*3;
@ -132,7 +156,7 @@ void ThirdOrder::command(int narg, char **arg)
int style = -1;
if (strcmp(arg[1],"regular") == 0) style = REGULAR;
else if (strcmp(arg[1],"eskm") == 0) style = BALLISTICO;
else if (strcmp(arg[1],"eskm") == 0) style = ESKM;
else error->all(FLERR,"Illegal Dynamical Matrix command");
// set option defaults
@ -143,15 +167,23 @@ void ThirdOrder::command(int narg, char **arg)
file_flag = 0;
file_opened = 0;
conversion = 1;
folded = 0;
// set Neigborlist attributes to NULL
ijnum = NULL;
neighbortags = NULL;
// read options from end of input line
if (style == REGULAR) options(narg-3,&arg[3]); //COME BACK
else if (style == BALLISTICO) options(narg-3,&arg[3]); //COME BACK
else if (comm->me == 0 && screen) fprintf(screen,"Illegal Dynamical Matrix command\n");
if (style == REGULAR) options(narg-3,&arg[3]);
else if (style == ESKM) options(narg-3,&arg[3]);
else error->all(FLERR,"Illegal Third Order command");
del = utils::numeric(FLERR, arg[2],false,lmp);
if (!folded) dynlenb = dynlen;
else dynlenb = (atom->natoms)*3;
if (atom->map_style == Atom::MAP_NONE)
error->all(FLERR,"third_order command requires an atom map, see atom_modify");
error->all(FLERR,"Third Order command requires an atom map, see atom_modify");
// move atoms by 3-vector or specified variable(s)
@ -163,7 +195,7 @@ void ThirdOrder::command(int narg, char **arg)
timer->barrier_stop();
}
if (style == BALLISTICO) {
if (style == ESKM) {
setup();
convert_units(update->unit_style);
conversion = conv_energy/conv_distance/conv_distance;
@ -183,13 +215,13 @@ void ThirdOrder::command(int narg, char **arg)
void ThirdOrder::options(int narg, char **arg)
{
if (narg < 0) error->all(FLERR,"Illegal third_order command");
if (narg < 0) error->all(FLERR,"Illegal Third Order command");
int iarg = 0;
const char *filename = "third_order.dat";
const char *filename = "Third Order.dat";
while (iarg < narg) {
if (strcmp(arg[iarg],"binary") == 0) {
if (iarg + 2 > narg) error->all(FLERR, "Illegal third_order command");
if (iarg + 2 > narg) error->all(FLERR, "Illegal Third Order command");
if (strcmp(arg[iarg+1],"gzip") == 0) {
compressed = 1;
} else {
@ -200,8 +232,15 @@ void ThirdOrder::options(int narg, char **arg)
if (iarg+2 > narg) error->all(FLERR, "Illegal third_order command");
filename = arg[iarg + 1];
file_flag = 1;
} else if (strcmp(arg[iarg],"fold") == 0) {
if (iarg+2 > narg) error->all(FLERR, "Illegal Third Order command");
if (strcmp(arg[iarg+1],"yes") == 0) {
folded = 1;
} else if (strcmp(arg[iarg+1],"no") == 0) {
folded = 0;
} else error->all(FLERR,"Illegal input for Third Order fold option");
iarg += 2;
} else error->all(FLERR,"Illegal third_order command");
} else error->all(FLERR,"Illegal Third Order command");
}
if (file_flag == 1 && me == 0) {
openfile(filename);
@ -231,11 +270,13 @@ void ThirdOrder::openfile(const char* filename)
}
if (!fp) error->one(FLERR,"Cannot open third_order file: {}", utils::getsyserror());
}
file_opened = 1;
}
/* ----------------------------------------------------------------------
create dynamical matrix
create third order tensor
------------------------------------------------------------------------- */
void ThirdOrder::calculateMatrix()
@ -247,28 +288,40 @@ void ThirdOrder::calculateMatrix()
bigint natoms = atom->natoms;
bigint *gm = groupmap;
double **f = atom->f;
int inum;
bigint j;
bigint *firstneigh;
double *dynmat = new double[3*dynlen];
double *fdynmat = new double[3*dynlen];
memset(&dynmat[0],0,dynlen*sizeof(double));
memset(&fdynmat[0],0,dynlen*sizeof(double));
double *dynmat = new double[dynlenb];
double *fdynmat = new double[dynlenb];
memset(&dynmat[0],0,dynlenb*sizeof(double));
memset(&fdynmat[0],0,dynlenb*sizeof(double));
getNeighbortags();
if (comm->me == 0 && screen) {
fprintf(screen,"Calculating Third Order ...\n");
fprintf(screen," Total # of atoms = " BIGINT_FORMAT "\n", natoms);
fprintf(screen," Atoms in group = " BIGINT_FORMAT "\n", gcount);
fprintf(screen," Total third order elements = "
BIGINT_FORMAT "\n", (dynlen*dynlen*dynlen) );
fprintf(screen, "Calculating Third Order ...\n");
fprintf(screen, " Total # of atoms = " BIGINT_FORMAT "\n", natoms);
fprintf(screen, " Atoms in group = " BIGINT_FORMAT "\n", gcount);
fprintf(screen, " Total third order elements = "
BIGINT_FORMAT "\n", (dynlen * dynlen * dynlen));
}
update->nsteps = 0;
int prog = 0;
for (bigint i=1; i<=natoms; i++) {
for (bigint i=1; i<=natoms; i++){
if (gm[i-1] < 0)
continue;
inum = ijnum[i-1];
firstneigh = neighbortags[i-1];
local_idx = atom->map(i);
for (int alpha=0; alpha<3; alpha++) {
for (bigint j=1; j<=natoms; j++) {
local_jdx = atom->map(j);
for (int beta=0; beta<3; beta++) {
for (int alpha=0; alpha<3; alpha++){
for (int jj=0; jj<inum; jj++){
j = firstneigh[jj];
if (gm[j] < 0 && !folded)
continue;
local_jdx = atom->map(j+1);
for (int beta=0; beta<3; beta++){
displace_atom(local_idx, alpha, 1);
displace_atom(local_jdx, beta, 1);
update_force();
@ -276,9 +329,13 @@ void ThirdOrder::calculateMatrix()
local_kdx = atom->map(k);
for (int gamma=0; gamma<3; gamma++) {
if (local_idx >= 0 && local_jdx >= 0 && local_kdx >= 0
&& gm[i-1] >= 0 && gm[j-1] >= 0 && gm[k-1] >= 0
&& ((gm[j] >= 0 && gm[k-1] >= 0) || folded)
&& local_kdx < nlocal) {
dynmat[gm[k-1]*3+gamma] += f[local_kdx][gamma];
if (folded) {
dynmat[(k-1)*3+gamma] += f[local_kdx][gamma];
} else {
dynmat[gm[k-1]*3+gamma] += f[local_kdx][gamma];
}
}
}
}
@ -288,9 +345,13 @@ void ThirdOrder::calculateMatrix()
local_kdx = atom->map(k);
for (int gamma=0; gamma<3; gamma++) {
if (local_idx >= 0 && local_jdx >= 0 && local_kdx >= 0
&& gm[i-1] >= 0 && gm[j-1] >= 0 && gm[k-1] >= 0
&& ((gm[j] >= 0 && gm[k-1] >= 0) || folded)
&& local_kdx < nlocal) {
dynmat[gm[k-1]*3+gamma] -= f[local_kdx][gamma];
if (folded) {
dynmat[(k-1)*3+gamma] -= f[local_kdx][gamma];
} else {
dynmat[gm[k-1]*3+gamma] -= f[local_kdx][gamma];
}
}
}
}
@ -302,9 +363,13 @@ void ThirdOrder::calculateMatrix()
local_kdx = atom->map(k);
for (int gamma=0; gamma<3; gamma++) {
if (local_idx >= 0 && local_jdx >= 0 && local_kdx >= 0
&& gm[i-1] >= 0 && gm[j-1] >= 0 && gm[k-1] >= 0
&& ((gm[j] >= 0 && gm[k-1] >= 0) || folded)
&& local_kdx < nlocal) {
dynmat[gm[k-1]*3+gamma] -= f[local_kdx][gamma];
if (folded) {
dynmat[(k-1)*3+gamma] -= f[local_kdx][gamma];
} else {
dynmat[gm[k-1]*3+gamma] -= f[local_kdx][gamma];
}
}
}
}
@ -314,24 +379,33 @@ void ThirdOrder::calculateMatrix()
local_kdx = atom->map(k);
for (int gamma=0; gamma<3; gamma++) {
if (local_idx >= 0 && local_jdx >= 0 && local_kdx >= 0
&& gm[i-1] >= 0 && gm[j-1] >= 0 && gm[k-1] >= 0
&& ((gm[j] >= 0 && gm[k-1] >= 0) || folded)
&& local_kdx < nlocal) {
dynmat[gm[k-1]*3+gamma] += f[local_kdx][gamma];
dynmat[gm[k-1]*3+gamma] /= (4 * del * del);
if (folded) {
dynmat[(k-1)*3+gamma] += f[local_kdx][gamma];
dynmat[(k-1)*3+gamma] /= (4 * del * del);
} else {
dynmat[gm[k-1]*3+gamma] += f[local_kdx][gamma];
dynmat[gm[k-1]*3+gamma] /= (4 * del * del);
}
}
}
}
displace_atom(local_jdx, beta, 1);
displace_atom(local_idx, alpha, 1);
MPI_Reduce(dynmat,fdynmat,3*dynlen,MPI_DOUBLE,MPI_SUM,0,world);
if (me == 0) {
writeMatrix(fdynmat, gm[i-1], alpha, gm[j-1], beta);
MPI_Reduce(dynmat,fdynmat,dynlenb,MPI_DOUBLE,MPI_SUM,0,world);
if (me == 0){
if (folded) {
writeMatrix(fdynmat, gm[i-1], alpha, j, beta);
} else {
writeMatrix(fdynmat, gm[i-1], alpha, gm[j], beta);
}
}
memset(&dynmat[0],0,dynlen*sizeof(double));
memset(&dynmat[0],0,dynlenb*sizeof(double));
}
}
}
if (comm->me == 0 && screen) {
if (me == 0 && screen) {
int p = 10 * gm[i-1] / gcount;
if (p > prog) {
prog = p;
@ -344,12 +418,12 @@ void ThirdOrder::calculateMatrix()
delete [] dynmat;
delete [] fdynmat;
if (screen && me ==0 )
if (screen && me ==0)
fprintf(screen,"Finished Calculating Third Order Tensor\n");
}
/* ----------------------------------------------------------------------
write dynamical matrix
write third order tensor
------------------------------------------------------------------------- */
void ThirdOrder::writeMatrix(double *dynmat, bigint i, int a, bigint j, int b)
@ -360,18 +434,34 @@ void ThirdOrder::writeMatrix(double *dynmat, bigint i, int a, bigint j, int b)
double norm;
if (!binaryflag && fp) {
clearerr(fp);
for (int k = 0; k < gcount; k++) {
norm = square(dynmat[k*3])+
square(dynmat[k*3+1])+
square(dynmat[k*3+2]);
if (norm > 1.0e-16)
fprintf(fp,
BIGINT_FORMAT " %d " BIGINT_FORMAT " %d " BIGINT_FORMAT
" %7.8f %7.8f %7.8f\n",
i+1, a + 1, j+1, b + 1, groupmap[k]+1,
dynmat[k*3] * conversion,
dynmat[k*3+1] * conversion,
dynmat[k*3+2] * conversion);
if (folded){
for (int k = 0; k < atom->natoms; k++){
norm = square(dynmat[k*3])+
square(dynmat[k*3+1])+
square(dynmat[k*3+2]);
if (norm > 1.0e-16)
fprintf(fp,
BIGINT_FORMAT " %d " BIGINT_FORMAT " %d " BIGINT_FORMAT
" %7.8f %7.8f %7.8f\n",
i+1, a + 1, j+1, b + 1, k+1,
dynmat[k*3] * conversion,
dynmat[k*3+1] * conversion,
dynmat[k*3+2] * conversion);
}
} else {
for (int k = 0; k < gcount; k++){
norm = square(dynmat[k*3])+
square(dynmat[k*3+1])+
square(dynmat[k*3+2]);
if (norm > 1.0e-16)
fprintf(fp,
BIGINT_FORMAT " %d " BIGINT_FORMAT " %d " BIGINT_FORMAT
" %7.8f %7.8f %7.8f\n",
i+1, a + 1, j+1, b + 1, groupmap[k]+1,
dynmat[k*3] * conversion,
dynmat[k*3+1] * conversion,
dynmat[k*3+2] * conversion);
}
}
} else if (binaryflag && fp) {
clearerr(fp);
@ -411,7 +501,18 @@ void ThirdOrder::displace_atom(int local_idx, int direction, int magnitude)
void ThirdOrder::update_force()
{
neighbor->ago = 0;
if ((modify->get_fix_by_id("package_intel")) ? true : false)
neighbor->decide();
force_clear();
int n_post_force = modify->n_post_force;
int n_pre_force = modify->n_pre_force;
int n_pre_reverse = modify->n_pre_reverse;
if (n_pre_force) {
modify->pre_force(vflag);
timer->stamp(Timer::MODIFY);
}
if (pair_compute_flag) {
force->pair->compute(eflag,vflag);
@ -428,10 +529,22 @@ void ThirdOrder::update_force()
force->kspace->compute(eflag,vflag);
timer->stamp(Timer::KSPACE);
}
if (n_pre_reverse) {
modify->pre_reverse(eflag,vflag);
timer->stamp(Timer::MODIFY);
}
if (force->newton) {
comm->reverse_comm();
timer->stamp(Timer::COMM);
}
// force modifications
if (n_post_force) {
modify->post_force(vflag);
timer->stamp(Timer::MODIFY);
}
++ update->nsteps;
}
@ -478,13 +591,13 @@ void ThirdOrder::convert_units(const char *style)
conv_distance = 1; // angstrom -> angstrom
} else if (strcmp(style,"si") == 0) {
if (comm->me) error->warning(FLERR,"Conversion Warning: Multiplication by Large Float");
if (me) error->warning(FLERR,"Conversion Warning: Multiplication by Large Float");
conv_energy = 6.022E22; // J -> 10 J/mol
conv_mass = 6.022E26; // kg -> g/mol
conv_distance = 1E-10; // meter -> angstrom
} else if (strcmp(style,"cgs") == 0) {
if (comm->me) error->warning(FLERR,"Conversion Warning: Multiplication by Large Float");
if (me) error->warning(FLERR,"Conversion Warning: Multiplication by Large Float");
conv_energy = 6.022E12; // Erg -> 10 J/mol
conv_mass = 6.022E23; // g -> g/mol
conv_distance = 1E-7; // centimeter -> angstrom
@ -495,13 +608,13 @@ void ThirdOrder::convert_units(const char *style)
conv_distance = 0.529177249; // bohr -> angstrom
} else if (strcmp(style,"micro") == 0) {
if (comm->me) error->warning(FLERR,"Conversion Warning: Untested Conversion");
if (me) error->warning(FLERR,"Conversion Warning: Untested Conversion");
conv_energy = 6.022E10; // picogram-micrometer^2/microsecond^2 -> 10 J/mol
conv_mass = 6.022E11; // pg -> g/mol
conv_distance = 1E-4; // micrometer -> angstrom
} else if (strcmp(style,"nano") == 0) {
if (comm->me) error->warning(FLERR,"Conversion Warning: Untested Conversion");
if (me) error->warning(FLERR,"Conversion Warning: Untested Conversion");
conv_energy = 6.022E4; // attogram-nanometer^2/nanosecond^2 -> 10 J/mol
conv_mass = 6.022E5; // ag -> g/mol
conv_distance = 0.1; // angstrom -> angstrom
@ -550,7 +663,7 @@ void ThirdOrder::create_groupmap()
for (int i=0; i<comm->nprocs; i++) {
recv[i] = 0;
}
recv[comm->me] = gid;
recv[me] = gid;
MPI_Allreduce(recv,displs,comm->nprocs,MPI_INT,MPI_SUM,world);
for (int i=0; i<comm->nprocs; i++) {
recv[i]=displs[i];
@ -579,3 +692,144 @@ void ThirdOrder::create_groupmap()
delete[] sub_groupmap;
delete[] temp_groupmap;
}
/* ---------------------------------------------------------------------- */
void ThirdOrder::getNeighbortags() {
// Create an extended neighbor list which is indexed by atom tag and yields atom tags
// groupmap[global atom index-1] = global atom indices (-1) of extended neighbors
bigint natoms = atom->natoms;
int *ilist,*jlist,*numneigh,**firstneigh;
bigint *Jlist,*klist;
int ii,jj,kk,inum,jnum,knum,sum;
int *temptags = (int*) malloc(natoms*sizeof(int));
int *ijnumproc = (int*) malloc(natoms*sizeof(int));
memory->create(ijnum, natoms, "thirdorder:ijnum");
bigint **firsttags;
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
memset(&ijnum[0],0,natoms*sizeof(int));
for (ii = 0; ii < inum; ii++) {
sum = 0;
memset(&temptags[0],0,natoms*sizeof(int));
jnum = numneigh[ii];
jlist = firstneigh[ii];
temptags[atom->tag[ilist[ii] & NEIGHMASK]-1] = 1;
for (jj = 0; jj < jnum; jj++) {
temptags[atom->tag[jlist[jj] & NEIGHMASK]-1] = 1;
}
for (bigint i=0; i<=natoms-1; i++) {
sum += temptags[i];
}
ijnum[atom->tag[ilist[ii] & NEIGHMASK]-1] = sum;
}
MPI_Allreduce(ijnum,ijnumproc,natoms,MPI_INT,MPI_SUM,world);
memset(&ijnum[0],0,natoms*sizeof(int));
sum = 0;
for (bigint i=0; i<=natoms-1; i++) {
sum += ijnumproc[i];
}
bigint nbytes = ((bigint) sizeof(bigint)) * sum;
bigint *data = (bigint *) memory->smalloc(nbytes, "thirdorder:firsttags");
bigint *datarecv = (bigint *) memory->smalloc(nbytes, "thirdorder:neighbortags");
nbytes = ((bigint) sizeof(bigint *)) * natoms;
firsttags = (bigint **) memory->smalloc(nbytes, "thirdorder:firsttags");
neighbortags = (bigint **) memory->smalloc(nbytes, "thirdorder:neighbortags");
memset(&data[0],0,sum*sizeof(bigint));
memset(&datarecv[0],0,sum*sizeof(bigint));
bigint n = 0;
for (bigint i = 0; i < natoms; i++) {
firsttags[i] = &data[n];
neighbortags[i] = &datarecv[n];
n += ijnumproc[i];
}
for (ii = 0; ii < inum; ii++) {
int m = 0;
memset(&temptags[0],0,natoms*sizeof(int));
jnum = numneigh[ii];
jlist = firstneigh[ii];
temptags[atom->tag[ilist[ii] & NEIGHMASK]-1] = 1;
for (jj = 0; jj < jnum; jj++) {
temptags[atom->tag[jlist[jj] & NEIGHMASK]-1] = 1;
}
for (int j=0; j < natoms; j++) {
if (temptags[j] == 1) {
neighbortags[atom->tag[ilist[ii] & NEIGHMASK]-1][m] = j;
m += 1;
}
}
}
MPI_Allreduce(datarecv,data,sum,MPI_LMP_BIGINT,MPI_SUM,world);
for (bigint i = 0; i < natoms; i++) {
ijnum[i] = 0;
sum = 0;
memset(&temptags[0],0,natoms*sizeof(int));
jnum = ijnumproc[i];
Jlist = firsttags[i];
temptags[i] = 1;
for (jj = 0; jj < jnum; jj++) {
temptags[Jlist[jj]] = 1;
klist = firsttags[Jlist[jj]];
knum = ijnumproc[Jlist[jj]];
for (kk = 0; kk < knum; kk++) {
temptags[klist[kk]] = 1;
}
}
for (bigint j=0; j<natoms; j++)
sum += temptags[j];
ijnum[i] = sum;
}
sum = 0;
for (bigint i=0; i<=natoms-1; i++) {
sum += ijnum[i];
}
free (neighbortags);
nbytes = ((bigint) sizeof(bigint)) * sum;
datarecv = (bigint *) memory->smalloc(nbytes, "thirdorder:firsttags");
nbytes = ((bigint) sizeof(bigint *)) * natoms;
neighbortags = (bigint **) memory->smalloc(nbytes, "thirdorder:neighbortags");
memset(&datarecv[0],0,sum*sizeof(bigint));
n = 0;
for (bigint i = 0; i < natoms; i++) {
neighbortags[i] = &datarecv[n];
n += ijnum[i];
}
for (bigint i = 0; i < natoms; i++) {
int m = 0;
memset(&temptags[0],0,natoms*sizeof(int));
jnum = ijnumproc[i];
Jlist = firsttags[i];
temptags[i] = 1;
for (int j = 0; j < jnum; j++) {
temptags[Jlist[j]] = 1;
klist = firsttags[Jlist[j]];
knum = ijnumproc[Jlist[j]];
for (kk = 0; kk < knum; kk++) {
temptags[klist[kk]] = 1;
}
}
for (bigint j=0; j < natoms; j++) {
if (temptags[j] == 1) {
neighbortags[i][m] = j;
m += 1;
}
}
}
free (firsttags);
free (ijnumproc);
free (temptags);
}

23
src/PHONON/third_order.h
View File

@ -23,8 +23,9 @@ class ThirdOrder : public Command {
void setup();
protected:
int eflag, vflag; // flags for energy/virial computation
int external_force_clear; // clear forces locally or externally
int eflag,vflag; // flags for energy/virial computation
int external_force_clear; // clear forces locally or externally
int triclinic; // 0 if domain is orthog, 1 if triclinic
int pairflag;
@ -34,25 +35,28 @@ class ThirdOrder : public Command {
int nvec; // local atomic dof = length of xvec
void update_force();
void force_clear();
virtual void openfile(const char *filename);
virtual void update_force();
virtual void force_clear();
virtual void openfile(const char* filename);
private:
protected:
void options(int, char **);
void create_groupmap();
void calculateMatrix();
void convert_units(const char *style);
void displace_atom(int local_idx, int direction, int magnitude);
void writeMatrix(double *, bigint, int, bigint, int);
void getNeighbortags();
double conversion;
double conv_energy;
double conv_distance;
double conv_mass;
double del;
int igroup, groupbit;
int igroup,groupbit;
bigint dynlen;
bigint dynlenb;
int scaleflag;
int me;
bigint gcount; // number of atoms in group
@ -62,6 +66,11 @@ class ThirdOrder : public Command {
int binaryflag; // 1 if dump file is written binary, 0 no
int file_opened; // 1 if openfile method has been called, 0 no
int file_flag; // 1 custom file name, 0 dynmat.dat
int folded; // 1 if system is folded, 0 no
class NeighList *list;
int *ijnum;
bigint **neighbortags;
FILE *fp;
};

16
src/input.cpp
View File

@ -784,9 +784,21 @@ int Input::execute_command()
if (flag) return 0;
// invoke commands added via style_command.h
// try suffixed version first
if (command_map->find(command) != command_map->end()) {
CommandCreator &command_creator = (*command_map)[command];
std::string mycmd = command;
if (lmp->suffix_enable) {
mycmd = command + std::string("/") + lmp->suffix;
if (command_map->find(mycmd) == command_map->end()) {
if (lmp->suffix2) {
mycmd = command + std::string("/") + lmp->suffix2;
if (command_map->find(mycmd) == command_map->end())
mycmd = command;
} else mycmd = command;
}
}
if (command_map->find(mycmd) != command_map->end()) {
CommandCreator &command_creator = (*command_map)[mycmd];
Command *cmd = command_creator(lmp);
cmd->command(narg,arg);
delete cmd;