diff --git a/doc/src/Commands_all.txt b/doc/src/Commands_all.txt index d0d7657c07..ecd21e42e7 100644 --- a/doc/src/Commands_all.txt +++ b/doc/src/Commands_all.txt @@ -124,6 +124,7 @@ An alphabetic list of all general LAMMPS commands. "thermo"_thermo.html, "thermo_modify"_thermo_modify.html, "thermo_style"_thermo_style.html, +"third_order"_third_order.html, "timer"_timer.html, "timestep"_timestep.html, "uncompute"_uncompute.html, diff --git a/doc/src/JPG/dynamical_matrix_dynmat.jpg b/doc/src/JPG/dynamical_matrix_dynmat.jpg new file mode 100644 index 0000000000..6c6adae72c Binary files /dev/null and b/doc/src/JPG/dynamical_matrix_dynmat.jpg differ diff --git a/doc/src/JPG/dynamical_matrix_force_constant.jpg b/doc/src/JPG/dynamical_matrix_force_constant.jpg new file mode 100644 index 0000000000..e1a0e50d14 Binary files /dev/null and b/doc/src/JPG/dynamical_matrix_force_constant.jpg differ diff --git a/doc/src/JPG/dynamical_matrix_phonons.jpg b/doc/src/JPG/dynamical_matrix_phonons.jpg new file mode 100644 index 0000000000..2a6d3a36d7 Binary files /dev/null and b/doc/src/JPG/dynamical_matrix_phonons.jpg differ diff --git a/doc/src/JPG/third_order_force_constant.png b/doc/src/JPG/third_order_force_constant.png new file mode 100644 index 0000000000..f6171ccf09 Binary files /dev/null and b/doc/src/JPG/third_order_force_constant.png differ diff --git a/doc/src/Packages_details.txt b/doc/src/Packages_details.txt index eead133add..fc30c045cf 100644 --- a/doc/src/Packages_details.txt +++ b/doc/src/Packages_details.txt @@ -1746,11 +1746,12 @@ USER-PHONON package :link(PKG-USER-PHONON),h4 A "fix phonon"_fix_phonon.html command that calculates dynamical matrices, which can then be used to compute phonon dispersion relations, directly from molecular dynamics simulations. -And a "dynamical_matrix" command to compute the dynamical matrix -from finite differences. +And a "dynamical_matrix"_dynamical_matrix.html as well as a +"third_order"_third_order.html command to compute the dynamical matrix +and third order tensor from finite differences. [Authors:] Ling-Ti Kong (Shanghai Jiao Tong University) for "fix phonon" -and Charlie Sievers (UC Davis) for "dynamical_matrix" +and Charlie Sievers (UC Davis) for "dynamical_matrix" and "third_order" [Supporting info:] @@ -1759,6 +1760,7 @@ src/USER-PHONON: filenames -> commands src/USER-PHONON/README "fix phonon"_fix_phonon.html "dynamical_matrix"_dynamical_matrix.html +"third_order"_third_order.html examples/USER/phonon :ul :line diff --git a/doc/src/commands_list.txt b/doc/src/commands_list.txt index a5c9b568ed..714aedefed 100644 --- a/doc/src/commands_list.txt +++ b/doc/src/commands_list.txt @@ -108,6 +108,7 @@ Commands :h1 thermo thermo_modify thermo_style + third_order timer timestep uncompute diff --git a/doc/src/dynamical_matrix.txt b/doc/src/dynamical_matrix.txt index f207297e9f..2d9b256694 100644 --- a/doc/src/dynamical_matrix.txt +++ b/doc/src/dynamical_matrix.txt @@ -30,14 +30,29 @@ dynamical_matrix 5 eskm 0.00000001 file dynamical.dat binary yes :pre [Description:] -Calculate the dynamical matrix of the selected group. +Calculate the dynamical matrix by finite difference of the selected group, + +:c,image(JPG/dynamical_matrix_dynmat.jpg) + +where D is the dynamical matrix and Phi is the force constant matrix defined by + +:c,image(JPG/dynamical_matrix_force_constant.jpg). + +The output for the dynamical matrix is printed three elements at a time. The +three elements are the three beta elements for a respective i/alpha/j combination. +Each line is printed in order of j increasing first, alpha second, and i last. + +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, where +frequencies, represented as omega, can be calculated from + +:c, image(Eqs/dynamical_matrix_phonons.jpg) [Restrictions:] -The command collects the entire dynamical matrix a single MPI rank, -so the memory requirements can be very significant for large systems. - -This command assumes a periodic system. +The command collects an array of nine times the number of atoms in a group +on every single MPI rank, so the memory requirements can be very significant +for large systems. This command is part of the USER-PHONON package. It is only enabled if LAMMPS was built with that package. See the "Build diff --git a/doc/src/lammps.book b/doc/src/lammps.book index 9868c8f299..eec7520fdc 100644 --- a/doc/src/lammps.book +++ b/doc/src/lammps.book @@ -217,6 +217,7 @@ temper_npt.html thermo.html thermo_modify.html thermo_style.html +third_order.html timer.html timestep.html uncompute.html diff --git a/doc/src/third_order.txt b/doc/src/third_order.txt new file mode 100644 index 0000000000..9636ec830e --- /dev/null +++ b/doc/src/third_order.txt @@ -0,0 +1,62 @@ +"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c + +:link(lws,http://lammps.sandia.gov) +:link(ld,Manual.html) +:link(lc,Commands_all.html) + +:line + +third_order command :h3 + +[Syntax:] + +third_order group-ID style delta args keyword value ... :pre + +group-ID = ID of group of atoms to displace :ulb,l +style = {regular} or {eskm} :l +delta = finite different displacement length (distance units) :l +one or more keyword/arg pairs may be appended :l + keyword = {file} or {binary} + {file} name = name of output file for the third order tensor + {binary} arg = {yes} or {no} or {gzip} :pre +:ule + +[Examples:] + +third_order 1 regular 0.000001 +third_order 1 eskm 0.000001 +third_order 3 regular 0.00004 file third_order.dat +third_order 5 eskm 0.00000001 file third_order.dat binary yes :pre + +[Description:] + +Calculate the third order force constant tensor by finite difference of the selected group, + +:c,image(JPG/third_order_force_constant.png)) + +where Phi is the third order force constant tensor. + +The output of the command is the tensor, three elements at a time. The +three elements correspond to the three gamma elements for a specific i/alpha/j/beta/k. +The initial five numbers are i, alpha, j, beta, and k respectively. + +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. + +[Restrictions:] + +The command collects a 9 times the number of atoms in the group on every single MPI rank, +so the memory requirements can be very significant for large systems. + +This command is part of the USER-PHONON package. It is only enabled if +LAMMPS was built with that package. See the "Build +package"_Build_package.html doc page for more info. + +[Related commands:] + +"fix phonon"_fix_phonon.html "dynamical_matrix"_dynamical_matrix.html + +[Default:] + +The default settings are file = "third_order.dat", binary = no diff --git a/doc/utils/sphinx-config/false_positives.txt b/doc/utils/sphinx-config/false_positives.txt index 13088139f9..bacce01a72 100644 --- a/doc/utils/sphinx-config/false_positives.txt +++ b/doc/utils/sphinx-config/false_positives.txt @@ -781,6 +781,7 @@ erotate Ertas ervel Espanol +eskm esu esub esw @@ -2129,6 +2130,7 @@ ph Phillpot phiphi phonon +phonons phophorous phosphide Phs diff --git a/examples/USER/phonon/dynamical_matrix_command/python/README.md b/examples/USER/phonon/dynamical_matrix_command/python/README.md new file mode 100755 index 0000000000..5b3c11febd --- /dev/null +++ b/examples/USER/phonon/dynamical_matrix_command/python/README.md @@ -0,0 +1,12 @@ +# LAMMPS LATTICE DYNAMICS COMMANDS + +## DYNAMICAL MATRIX CALCULATOR + +This directory contains the ingredients to calculate a dynamical matrix with python. + +Example: +``` +python dynmat.py +``` + +## Requires: MANYBODY and MOLECULE packages and the Python Library Interface diff --git a/examples/USER/phonon/dynamical_matrix_command/python/dynmat.py b/examples/USER/phonon/dynamical_matrix_command/python/dynmat.py new file mode 100644 index 0000000000..2a3a0b5a2f --- /dev/null +++ b/examples/USER/phonon/dynamical_matrix_command/python/dynmat.py @@ -0,0 +1,42 @@ +"""Made by Charlie Sievers Ph.D. Candidate, UC Davis, Donadio Lab 2019""" +# from mpi4py import MPI +from lammps import lammps +import numpy as np + +# comm = MPI.COMM_WORLD +# rank = comm.Get_rank() + +""" LAMMPS VARIABLES """ + +# data files +infile = "silicon_input_file.lmp" +ff_file = "ff-silicon.lmp" + +# full output useful for testing +lmp = lammps() + +# reduced output useful reducing IO for production runs +# lmp = lammps(cmdargs=["-screen", "none", "-log", "none"]) + +# lammps commands +lmp.command("atom_style full") +lmp.command("units metal") +lmp.command("processors * * *") +lmp.command("neighbor 1 bin") +lmp.command("boundary p p p") + +# read data and force field file +lmp.command("read_data {}".format(infile)) +lmp.file("{}".format(ff_file)) + +lmp.command("dynamical_matrix all eskm 0.000001 file dynmat.dat") + +dynmat = np.loadtxt("dynmat.dat") +dynlen = int(3*np.sqrt(len(dynmat)/3)) +dynmat = dynmat.reshape((dynlen, dynlen)) + +eigvals, eigvecs = np.linalg.eig(dynmat) + +# frequencies in THz +omegas = np.sqrt(np.abs(eigvals)) +print(omegas) diff --git a/src/USER-PHONON/Install.sh b/src/USER-PHONON/Install.sh index 26104b45cf..a73f529cfa 100755 --- a/src/USER-PHONON/Install.sh +++ b/src/USER-PHONON/Install.sh @@ -42,3 +42,5 @@ action fix_phonon.cpp fft3d_wrap.h action fix_phonon.h fft3d_wrap.h action dynamical_matrix.cpp action dynamical_matrix.h +action third_order.cpp +action third_order.h diff --git a/src/USER-PHONON/README b/src/USER-PHONON/README index b554eacd5e..d5ed666c0c 100644 --- a/src/USER-PHONON/README +++ b/src/USER-PHONON/README @@ -3,11 +3,11 @@ matrices from finite temperature MD simulations, which can then be used to compute phonon dispersion relations, directly from molecular dynamics simulations. -It also contains a command to compute the dynamical matrix at -pre-optimized positions through finite differences. +It also contains commands to compute the dynamical matrix and third +order tensor at pre-optimized positions through finite differences. See the doc page for the fix phonon command or the dynamical_matrix -command for detailed usage instructions. +or third_order commands for detailed usage instructions. Use of this package requires building LAMMPS with FFT suppport, as described in doc/Section_start.html. diff --git a/src/USER-PHONON/dynamical_matrix.cpp b/src/USER-PHONON/dynamical_matrix.cpp index 559ef4c36f..fe266fba76 100644 --- a/src/USER-PHONON/dynamical_matrix.cpp +++ b/src/USER-PHONON/dynamical_matrix.cpp @@ -66,11 +66,6 @@ void DynamicalMatrix::setup() domain->image_check(); domain->box_too_small_check(); neighbor->build(1); - neighbor->ncalls = 0; - neighbor->every = 2; // build every this many steps - neighbor->delay = 1; - neighbor->ago = 0; - neighbor->ndanger = 0; // compute all forces external_force_clear = 0; @@ -273,7 +268,7 @@ void DynamicalMatrix::calculateMatrix() local_idx = atom->map(i); if (gm[i-1] < 0) continue; - for (bigint alpha=0; alpha<3; alpha++){ + for (int alpha=0; alpha<3; alpha++){ displace_atom(local_idx, alpha, 1); update_force(); for (bigint j=1; j<=natoms; j++){ @@ -291,7 +286,7 @@ void DynamicalMatrix::calculateMatrix() local_jdx = atom->map(j); if (local_idx >= 0 && local_jdx >= 0 && local_jdx < nlocal && gm[j-1] >= 0){ - for (bigint beta=0; beta<3; beta++){ + for (int beta=0; beta<3; beta++){ if (atom->rmass_flag == 1) imass = sqrt(m[local_idx] * m[local_jdx]); else diff --git a/src/USER-PHONON/third_order.cpp b/src/USER-PHONON/third_order.cpp new file mode 100644 index 0000000000..7764287337 --- /dev/null +++ b/src/USER-PHONON/third_order.cpp @@ -0,0 +1,575 @@ +// +// Created by charlie sievers on 7/5/18. +// + +#include "third_order.h" +#include +#include +#include +#include "atom.h" +#include "domain.h" +#include "comm.h" +#include "error.h" +#include "group.h" +#include "force.h" +#include "memory.h" +#include "bond.h" +#include "angle.h" +#include "dihedral.h" +#include "improper.h" +#include "kspace.h" +#include "update.h" +#include "neighbor.h" +#include "pair.h" +#include "timer.h" +#include "finish.h" +#include "math_special.h" +#include +#include + +using namespace LAMMPS_NS; +using namespace MathSpecial; +enum{REGULAR,BALLISTICO}; + +/* ---------------------------------------------------------------------- */ + +ThirdOrder::ThirdOrder(LAMMPS *lmp) : Pointers(lmp), fp(NULL) +{ + external_force_clear = 1; +} + +/* ---------------------------------------------------------------------- */ + +ThirdOrder::~ThirdOrder() +{ + if (fp && me == 0) fclose(fp); + fp = NULL; + memory->destroy(groupmap); +} + +/* ---------------------------------------------------------------------- + setup without output or one-time post-init setup + flag = 0 = just force calculation + flag = 1 = reneighbor and force calculation +------------------------------------------------------------------------- */ + +void ThirdOrder::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); + + // compute all forces + external_force_clear = 0; + eflag=0; + vflag=0; + update_force(); + + if (gcount == atom->natoms) + for (bigint i=0; inatoms; i++) + groupmap[i] = i; + else + create_groupmap(); +} + +/* ---------------------------------------------------------------------- */ + +void ThirdOrder::command(int narg, char **arg) +{ + MPI_Comm_rank(world,&me); + + if (domain->box_exist == 0) + error->all(FLERR,"third_order command before simulation box is defined"); + if (narg < 2) error->all(FLERR,"Illegal third_order command"); + + lmp->init(); + + // orthogonal vs triclinic simulation box + + 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; + + // 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; + + int style = -1; + if (strcmp(arg[1],"regular") == 0) style = REGULAR; + else if (strcmp(arg[1],"eskm") == 0) style = BALLISTICO; + else error->all(FLERR,"Illegal Dynamical Matrix command"); + + // set option defaults + + binaryflag = 0; + scaleflag = 0; + compressed = 0; + file_flag = 0; + file_opened = 0; + conversion = 1; + + // 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"); + del = force->numeric(FLERR, arg[2]); + + if (atom->map_style == 0) + error->all(FLERR,"third_order command requires an atom map, see atom_modify"); + + // move atoms by 3-vector or specified variable(s) + + if (style == REGULAR) { + setup(); + timer->init(); + timer->barrier_start(); + calculateMatrix(); + timer->barrier_stop(); + } + + if (style == BALLISTICO) { + setup(); + convert_units(update->unit_style); + conversion = conv_energy/conv_distance/conv_distance; + timer->init(); + timer->barrier_start(); + calculateMatrix(); + timer->barrier_stop(); + } + + Finish finish(lmp); + finish.end(1); +} + +/* ---------------------------------------------------------------------- + parse optional parameters +------------------------------------------------------------------------- */ + +void ThirdOrder::options(int narg, char **arg) +{ + if (narg < 0) error->all(FLERR,"Illegal third_order command"); + int iarg = 0; + const char *filename = "third_order.dat"; + std::stringstream fss; + + while (iarg < narg) { + if (strcmp(arg[iarg],"file") == 0) { + if (iarg+2 > narg) error->all(FLERR, "Illegal third_order command"); + fss << arg[iarg + 1]; + filename = fss.str().c_str(); + file_flag = 1; + iarg += 2; + } else if (strcmp(arg[iarg],"binary") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal third_order command"); + if (strcmp(arg[iarg+1],"gzip") == 0) { + compressed = 1; + } else if (strcmp(arg[iarg+1],"yes") == 0) { + binaryflag = 1; + } + iarg += 2; + } else error->all(FLERR,"Illegal third_order command"); + } + if (file_flag == 1 and me == 0) { + openfile(filename); + } +} + +/* ---------------------------------------------------------------------- + generic opening of a file + ASCII or binary or gzipped + some derived classes override this function +------------------------------------------------------------------------- */ + +void ThirdOrder::openfile(const char* filename) +{ + // if file already opened, return + if (file_opened) return; + + if (compressed) { +#ifdef LAMMPS_GZIP + char gzip[128]; + sprintf(gzip,"gzip -6 > %s",filename); +#ifdef _WIN32 + fp = _popen(gzip,"wb"); +#else + fp = popen(gzip,"w"); +#endif +#else + error->one(FLERR,"Cannot open gzipped file"); +#endif + } else if (binaryflag) { + fp = fopen(filename,"wb"); + } else { + fp = fopen(filename,"w"); + } + + if (fp == NULL) error->one(FLERR,"Cannot open dump file"); + + file_opened = 1; +} + +/* ---------------------------------------------------------------------- + create dynamical matrix +------------------------------------------------------------------------- */ + +void ThirdOrder::calculateMatrix() +{ + int local_idx; // local index + int local_jdx; // second local index + int local_kdx; // third local index + int nlocal = atom->nlocal; + bigint natoms = atom->natoms; + bigint *gm = groupmap; + double **f = atom->f; + + 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)); + + 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) ); + } + + update->nsteps = 0; + int prog = 0; + for (bigint i=1; i<=natoms; i++){ + 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++){ + displace_atom(local_idx, alpha, 1); + displace_atom(local_jdx, beta, 1); + update_force(); + for (bigint k=1; k<=natoms; k++){ + 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 + && local_kdx < nlocal) { + dynmat[gm[k-1]*3+gamma] += f[local_kdx][gamma]; + } + } + } + displace_atom(local_jdx, beta, -2); + update_force(); + for (bigint k=1; k<=natoms; k++){ + 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 + && local_kdx < nlocal) { + dynmat[gm[k-1]*3+gamma] -= f[local_kdx][gamma]; + } + } + } + displace_atom(local_jdx, beta, 1); + displace_atom(local_idx,alpha,-2); + displace_atom(local_jdx, beta, 1); + update_force(); + for (bigint k=1; k<=natoms; k++){ + 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 + && local_kdx < nlocal) { + dynmat[gm[k-1]*3+gamma] -= f[local_kdx][gamma]; + } + } + } + displace_atom(local_jdx, beta, -2); + update_force(); + for (bigint k=1; k<=natoms; k++){ + 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 + && local_kdx < nlocal) { + 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); + } + memset(&dynmat[0],0,dynlen*sizeof(double)); + } + } + } + if (comm->me == 0 && screen) { + int p = 10 * gm[i-1] / gcount; + if (p > prog) { + prog = p; + fprintf(screen," %d%%",p*10); + fflush(screen); + } + } + } + + delete [] dynmat; + delete [] fdynmat; + + if (screen && me ==0 ) + fprintf(screen,"Finished Calculating Third Order Tensor\n"); +} + +/* ---------------------------------------------------------------------- + write dynamical matrix +------------------------------------------------------------------------- */ + +void ThirdOrder::writeMatrix(double *dynmat, bigint i, int a, bigint j, int b) +{ + if (me != 0) + return; + + 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); + } + } else if (binaryflag && fp){ + clearerr(fp); + fwrite(&dynmat[0], sizeof(double), dynlen, fp); + } + if (ferror(fp)) error->one(FLERR,"Error writing to file"); + +} + +/* ---------------------------------------------------------------------- + Displace atoms + ---------------------------------------------------------------------- */ + +void ThirdOrder::displace_atom(int local_idx, int direction, int magnitude) +{ + if (local_idx < 0) return; + + double **x = atom->x; + int *sametag = atom->sametag; + int j = local_idx; + + x[local_idx][direction] += del*magnitude; + + while (sametag[j] >= 0){ + j = sametag[j]; + x[j][direction] += del*magnitude; + } +} + +/* ---------------------------------------------------------------------- + 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 ThirdOrder::update_force() +{ + force_clear(); + + if (pair_compute_flag) { + force->pair->compute(eflag,vflag); + timer->stamp(Timer::PAIR); + } + if (atom->molecular) { + 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); + } + ++ update->nsteps; +} + +/* ---------------------------------------------------------------------- + clear force on own & ghost atoms + clear other arrays as needed +------------------------------------------------------------------------- */ + +void ThirdOrder::force_clear() +{ + if (external_force_clear) return; + + // 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; + + if (nbytes) { + memset(&atom->f[0][0],0,3*nbytes); + } +} + +/* ---------------------------------------------------------------------- */ + +void ThirdOrder::convert_units(const char *style) +{ + // physical constants from: + // http://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 + + } 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,"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,"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,"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,"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 error->all(FLERR,"Units Type Conversion Not Found"); + +} + +/* ---------------------------------------------------------------------- */ + +void ThirdOrder::create_groupmap() +{ + //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]; + + //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; + } + } + + //populate arrays for Allgatherv + for (int i=0; inprocs; i++){ + recv[i] = 0; + } + recv[comm->me] = gid; + MPI_Allreduce(recv,displs,comm->nprocs,MPI_INT,MPI_SUM,world); + for (int i=0; inprocs; 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); + + //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; +} diff --git a/src/USER-PHONON/third_order.h b/src/USER-PHONON/third_order.h new file mode 100644 index 0000000000..83062b6b1f --- /dev/null +++ b/src/USER-PHONON/third_order.h @@ -0,0 +1,76 @@ +// +// Created by charlie sievers on 7/5/18. +// + + +#ifdef COMMAND_CLASS + +CommandStyle(third_order,ThirdOrder) + +#else + +#ifndef LMP_THIRD_ORDER_H +#define LMP_THIRD_ORDER_H + +#include "pointers.h" + +namespace LAMMPS_NS { + + class ThirdOrder : protected Pointers { + public: + ThirdOrder(class LAMMPS *); + virtual ~ThirdOrder(); + void command(int, char **); + void setup(); + + protected: + 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; + + int pair_compute_flag; // 0 if pair->compute is skipped + int kspace_compute_flag; // 0 if kspace->compute is skipped + + int nvec; // local atomic dof = length of xvec + + void update_force(); + void force_clear(); + virtual void openfile(const char* filename); + + + private: + 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); + + double conversion; + double conv_energy; + double conv_distance; + double conv_mass; + double del; + int igroup,groupbit; + bigint dynlen; + int scaleflag; + int me; + bigint gcount; // number of atoms in group + bigint *groupmap; + + int compressed; // 1 if dump file is written compressed, 0 no + 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 + + FILE *fp; + }; +} + + +#endif //LMP_THIRD_ORDER_H +#endif +