diff --git a/doc/src/compute_born.rst b/doc/src/compute_born.rst new file mode 100644 index 0000000000..c567a96d62 --- /dev/null +++ b/doc/src/compute_born.rst @@ -0,0 +1,130 @@ +.. index:: compute born + +compute born command +==================== + +Syntax +"""""" + +.. parsed-literal:: + + compute ID group-ID born + +* ID, group-ID are documented in :doc:`compute ` command +* born = style name of this compute command + +Examples +"""""""" + +.. code-block:: LAMMPS + + compute 1 all born + +Description +""""""""""" + +Define a compute that calculates +:math:`\frac{\partial{}^2U}{\partial\varepsilon_{i}\partial\varepsilon_{j}}` the +second derivatives of the potential energy :math:`U` with regard to strain +tensor :math:`\varepsilon` elements. These values are related to: + +.. math:: + + C^{B}_{i,j}=\frac{1}{V}\frac{\partial{}^2U}{\partial{}\varepsilon_{i}\partial\varepsilon_{j}} + +also called the Born term of elastic contants in the stress-stress fluctuation +formalism. This quantity can be used to compute the elastic constant tensor. +Using the symmetric Voigt notation, the elastic constant tensor can be written +as a 6x6 symmetric matrix: + +.. math:: + + C_{i,j} = \langle{}C^{B}_{i,j}\rangle + + \frac{V}{k_{B}T}\left(\langle\sigma_{i}\sigma_{j}\rangle\right. + \left.- \langle\sigma_{i}\rangle\langle\sigma_{j}\rangle\right) + + \frac{Nk_{B}T}{V} + \left(\delta_{i,j}+(\delta_{1,i}+\delta_{2,i}+\delta_{3,i})\right. + \left.*(\delta_{1,j}+\delta_{2,j}+\delta_{3,j})\right) + +In the above expression, :math:`\sigma` stands for the virial stress +tensor, :math:`\delta` is the Kronecker delta and the usual notation apply for +the number of particle, the temperature and volume respectively :math:`N`, +:math:`T` and :math:`V`. :math `k_{B}` is the Boltzmann constant. + +The Born term is a symmetric 6x6 matrix by construction and as such can be +expressed as 21 independent terms. The terms are ordered corresponding to the +following matrix element: + +.. math:: + + \matrix{ + C_{1} & C_{7} & C_{8} & C_{9} & C_{10} & C_{11} \\ + C_{7} & C_{2} & C_{12} & C_{13} & C_{14} & C_{15} \\ + \vdots & C_{12} & C_{3} & C_{16} & C_{17} & C_{18} \\ + \vdots & C_{13} & C_{16} & C_{4} & C_{19} & C_{20} \\ + \vdots & \vdots & \vdots & C_{19} & C_{5} & C_{21} \\ + \vdots & \vdots & \vdots & \vdots & C_{21} & C_{6} + } + +in this matrix the indices of :math:`C_{k}` value are the corresponding index +:math:`k` in the compute output. Each term comes from the sum of every +interactions derivatives in the system as explained in :ref:`(VanWorkum) +` or :ref:`(Voyiatzis) `. + +The output can be accessed using usual Lammps routines: + +.. code-block:: LAMMPS + + compute 1 all born + compute 2 all pressure NULL virial + variable S1 equal -c_2[1] + variable S2 equal -c_2[2] + variable S3 equal -c_2[3] + variable S4 equal -c_2[4] + variable S5 equal -c_2[5] + variable S6 equal -c_2[6] + fix 1 all ave/time 1 1 1 v_S1 v_S2 v_S3 v_S4 v_S5 v_S6 c_1[*] file born.out + +In this example, the file *born.out* will contain the information needed to +compute the first and second terms of the elastic constant matrix in a post +processing procedure. The other required quantities can be accessed using any +other *LAMMPS* usual method. + +NOTE: In the above :math:`C_{i,j}` computation, the term involving the virial +stress tensor :math:`\sigma` is the covariance between each elements. In a +solid the virial stress can have large variations between timesteps and average +values can be slow to converge. This term is better computed using +instantaneous values. + +**Output info:** + +This compute calculates a global array with the number of rows=21. +The values are ordered as explained above. These values can be used +by any command that uses a global values from a compute as input. See +the :doc:`Howto output ` doc page for an overview of +LAMMPS output options. + +The array values calculated by this compute are all "extensive". + +Restrictions +"""""""""""" + +The Born term can be decomposed as a product of two terms. The first one +is a general term which depends on the configuration. The second one is +specific to every interaction composing your forcefield (non-bonded, +bonds, angle...). Currently not all interaction implement the born +method giving first and second order derivatives and an error will +be raised if you try to use this compute with such interactions. + +Default +""""""" + +none + +.. _VanWorkum: + +K.Van Workum et al. J. Chem. Phys. 125 144506 (2006) + +.. _Voyiatzis: + +E.Voyiatzis, Computer Physics Communications 184(2013)27-33 diff --git a/src/angle.cpp b/src/angle.cpp index d822418141..969ae898bb 100644 --- a/src/angle.cpp +++ b/src/angle.cpp @@ -35,6 +35,7 @@ Angle::Angle(LAMMPS *lmp) : Pointers(lmp) allocated = 0; suffix_flag = Suffix::NONE; + born_enable = 0; maxeatom = maxvatom = maxcvatom = 0; eatom = nullptr; diff --git a/src/angle.h b/src/angle.h index 3a6521003e..1c7f4e6fc3 100644 --- a/src/angle.h +++ b/src/angle.h @@ -26,6 +26,7 @@ class Angle : protected Pointers { int allocated; int *setflag; int writedata; // 1 if writes coeffs to data file + int born_enable; double energy; // accumulated energies double virial[6]; // accumulated virial: xx,yy,zz,xy,xz,yz double *eatom, **vatom; // accumulated per-atom energy/virial @@ -56,6 +57,11 @@ class Angle : protected Pointers { virtual void read_restart_settings(FILE *){}; virtual void write_data(FILE *) {} virtual double single(int, int, int, int) = 0; + virtual void born(int/*atype*/, int/*at1*/, int/*at2*/, int/*at3*/, double& du, double& du2) + { + du = 0.0; + du2 = 0.0; + } virtual double memory_usage(); protected: diff --git a/src/bond.cpp b/src/bond.cpp index fb313f0170..6d037add53 100644 --- a/src/bond.cpp +++ b/src/bond.cpp @@ -41,6 +41,7 @@ Bond::Bond(LAMMPS *lmp) : Pointers(lmp) allocated = 0; suffix_flag = Suffix::NONE; + born_enable = 0; maxeatom = maxvatom = 0; eatom = nullptr; diff --git a/src/bond.h b/src/bond.h index b2b2008b9e..402a4e644d 100644 --- a/src/bond.h +++ b/src/bond.h @@ -30,6 +30,7 @@ class Bond : protected Pointers { double virial[6]; // accumulated virial: xx,yy,zz,xy,xz,yz double *eatom, **vatom; // accumulated per-atom energy/virial + int born_enable; int reinitflag; // 1 if compatible with fix adapt and alike // KOKKOS host/device flag and data masks @@ -56,6 +57,12 @@ class Bond : protected Pointers { virtual void *extract(const char *, int &) { return nullptr; } virtual void reinit(); + virtual void born(int/*btype*/, double/*rsq*/, int/*at1*/, int/*at2*/, double& du, double& du2) + { + du = 0.0; + du2 = 0.0; + } + void write_file(int, char **); protected: diff --git a/src/compute_born.cpp b/src/compute_born.cpp new file mode 100644 index 0000000000..76059512a6 --- /dev/null +++ b/src/compute_born.cpp @@ -0,0 +1,806 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. + ------------------------------------------------------------------------- */ + +/*------------------------------------------------------------------------ + Contributing Authors : Germain Clavier (UCA) + --------------------------------------------------------------------------*/ + +#include "compute_born.h" +#include +#include +#include + +#include "atom.h" +#include "atom_vec.h" +#include "update.h" +#include "domain.h" +#include "neighbor.h" +#include "force.h" +#include "pair.h" +#include "bond.h" +#include "angle.h" +#include "dihedral.h" +#include "improper.h" +#include "molecule.h" +#include "neigh_request.h" +#include "neigh_list.h" +#include "error.h" +#include "memory.h" + +using namespace LAMMPS_NS; + +#define BIG 1000000000 + + +/* ---------------------------------------------------------------------- */ + +ComputeBorn::ComputeBorn(LAMMPS *lmp, int narg, char **arg) : + Compute(lmp, narg, arg) +{ + MPI_Comm_rank(world,&me); + + // For now the matrix can be computed as a 21 element vector + + nvalues = 21; + + // Error check + + // Initialize some variables + + values_local = values_global = vector = NULL; + + // this fix produces a global vector + + memory->create(vector,nvalues,"born:vector"); + memory->create(values_local,nvalues,"born:values_local"); + memory->create(values_global,nvalues,"born:values_global"); + size_vector = nvalues; + + vector_flag = 1; + extvector = 0; + +} + +/* ---------------------------------------------------------------------- */ + +ComputeBorn::~ComputeBorn() +{ + + // delete [] which; + + memory->destroy(values_local); + memory->destroy(values_global); + memory->destroy(vector); +} + +/* ---------------------------------------------------------------------- */ + +void ComputeBorn::init() +{ + + // Timestep Value + dt = update->dt; + + pairflag = 0; + bondflag = 0; + angleflag = 0; + dihedflag = 0; + impflag = 0; + + // Error check + + // For now this compute requires at least a pair style with pair_born method + // implemented + + if (force->pair == NULL) + error->all(FLERR,"No pair style is defined for compute born"); + if (force->pair->born_enable == 0) { + error->all(FLERR,"Pair style does not support compute born"); + pairflag = 0; + } else { + pairflag = 1; + } + + // Throws an error for now? + if (force->bond != NULL) { + if (force->bond->born_enable == 0) { + error->warning(FLERR, "Bond style does not support compute born"); + bondflag = 0; + } else { + bondflag = 1; + } + } + + if (force->angle != NULL) { + if (force->angle->born_enable == 0) { + error->warning(FLERR, "Angle style does not support compute born"); + angleflag = 0; + } else { + angleflag = 1; + } + } + + if (force->dihedral != NULL) { + if (force->dihedral->born_enable == 0) { + error->warning(FLERR, "Dihedral style does not support compute born"); + dihedflag = 0; + } else { + dihedflag = 1; + } + } + + if (force->improper != NULL) { + if (force->improper->born_enable == 0) { + error->warning(FLERR, "Improper style does not support compute born"); + impflag = 0; + } else { + impflag = 1; + } + } + + // need an occasional half neighbor list + int irequest = neighbor->request((void *) this); + neighbor->requests[irequest]->pair = 0; + neighbor->requests[irequest]->compute = 1; + neighbor->requests[irequest]->occasional = 1; +} + +/* ---------------------------------------------------------------------- */ + +void ComputeBorn::init_list(int /* id */, NeighList *ptr) +{ + list = ptr; +} + +/* ---------------------------------------------------------------------- + compute output vector + ------------------------------------------------------------------------- */ + +void ComputeBorn::compute_vector() +{ + invoked_array = update->ntimestep; + + // zero out arrays for one sample + + int i; + for (i = 0; i < nvalues; i++) values_local[i] = 0.0; + + // Compute Born contribution + if (pairflag) compute_pairs(); + + // For now these functions are commented + if (bondflag) compute_bonds(); + + if (angleflag) compute_angles(); + + if (dihedflag) compute_dihedrals(); + + // Even if stated in Voyatzis-2012, improper and dihedrals + // are not exactly the same in lammps. Atoms order can depend + // on the forcefield/improper interaction used. As such, + // writing a general routine to compute improper contribution + // might be more tricky than expected. + // if (impflag) compute_impropers(); + + // sum Born contributions over all procs + MPI_Allreduce(values_local,values_global,nvalues, + MPI_DOUBLE,MPI_SUM,world); + + int m; + for (m=0; mtype; + int *mask = atom->mask; + int nlocal = atom->nlocal; + double *special_coul = force->special_coul; + double *special_lj = force->special_lj; + int newton_pair = force->newton_pair; + + // invoke half neighbor list (will copy or build if necessary) + neighbor->build_one(list); + + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + + // loop over neighbors of my atoms + + Pair *pair = force->pair; + double **cutsq = force->pair->cutsq; + + // Declares born values + + int a,b,c,d; + double xi[3]; + double fi[3]; + double xj[3]; + double rij[3]; + double pair_pref; + double r2inv; + + + m = 0; + while (mx[i][0]; + xi[1] = atom->x[i][1]; + xi[2] = atom->x[i][2]; + itype = type[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + factor_lj = special_lj[sbmask(j)]; + factor_coul = special_coul[sbmask(j)]; + j &= NEIGHMASK; + + if (!(mask[j] & groupbit)) continue; + + xj[0] = atom->x[j][0]; + xj[1] = atom->x[j][1]; + xj[2] = atom->x[j][2]; + delx = xi[0] - xj[0]; + dely = xi[1] - xj[1]; + delz = xi[2] - xj[2]; + rij[0] = xj[0]-xi[0]; + rij[1] = xj[1]-xi[1]; + rij[2] = xj[2]-xi[2]; + rsq = delx*delx + dely*dely + delz*delz; + r2inv = 1.0/rsq; + jtype = type[j]; + + if (rsq >= cutsq[itype][jtype]) continue; + + if (newton_pair || j < nlocal) { + // Add contribution to Born tensor + + pair->born(i,j,itype,jtype,rsq,factor_coul,factor_lj,dupair,du2pair); + pair_pref = du2pair - dupair*rinv; + + // See albemunu in compute_born.h for indices order. + a = 0; + b = 0; + c = 0; + d = 0; + for (i = 0; i<21; i++) { + a = albemunu[i][0]; + b = albemunu[i][1]; + c = albemunu[i][2]; + d = albemunu[i][3]; + values_local[m+i] += pair_pref*rij[a]*rij[b]*rij[c]*rij[d]*r2inv; + } + } + } + } + m += 21; + } +} + +/* ---------------------------------------------------------------------- + count bonds and compute bond info on this proc + only count bond once if newton_bond is off + all atoms in interaction must be in group + all atoms in interaction must be known to proc + if bond is deleted or turned off (type <= 0) + do not count or count contribution + COMMENTED FOR NOW +---------------------------------------------------------------------- */ +void ComputeBorn::compute_bonds() +{ +/* ---------------------------------------------------------------------- + int i,m,n,nb,atom1,atom2,imol,iatom,btype,ivar; + tagint tagprev; + double dx,dy,dz,rsq; + + double **x = atom->x; + double **v = atom->v; + int *type = atom->type; + tagint *tag = atom->tag; + int *num_bond = atom->num_bond; + tagint **bond_atom = atom->bond_atom; + int **bond_type = atom->bond_type; + int *mask = atom->mask; + + int *molindex = atom->molindex; + int *molatom = atom->molatom; + Molecule **onemols = atom->avec->onemols; + + int nlocal = atom->nlocal; + int newton_bond = force->newton_bond; + int molecular = atom->molecular; + + Bond *bond = force->bond; + + int a,b,c,d; + double rij[3]; + double rinv, r2inv; + double pair_pref, dupair, du2pair; + + // loop over all atoms and their bonds + + m = 0; + while (mnum_bond[iatom]; + } + + for (i = 0; i < nb; i++) { + if (molecular == 1) { + btype = bond_type[atom1][i]; + atom2 = atom->map(bond_atom[atom1][i]); + } else { + tagprev = tag[atom1] - iatom - 1; + btype = onemols[imol]->bond_type[iatom][i]; + atom2 = atom->map(onemols[imol]->bond_atom[iatom][i]+tagprev); + } + + if (atom2 < 0 || !(mask[atom2] & groupbit)) continue; + if (newton_bond == 0 && tag[atom1] > tag[atom2]) continue; + if (btype <= 0) continue; + + dx = x[atom2][0] - x[atom1][0]; + dy = x[atom2][1] - x[atom1][1]; + dz = x[atom2][2] - x[atom1][2]; + domain->minimum_image(dx,dy,dz); + rsq = dx*dx + dy*dy + dz*dz; + rij[0] = dx; + rij[1] = dy; + rij[2] = dz; + r2inv = 1.0/rsq; + rinv = sqrt(r2inv); + + pair_pref = 0.0; + dupair = 0.0; + du2pair = 0.0; + bond->born(btype,rsq,atom1,atom2,dupair,du2pair); + pair_pref = du2pair - dupair*rinv; + + a = 0; + b = 0; + c = 0; + d = 0; + for (i = 0; i<21; i++) { + a = albemunu[i][0]; + b = albemunu[i][1]; + c = albemunu[i][2]; + d = albemunu[i][3]; + values_local[m+i] += pair_pref*rij[a]*rij[b]*rij[c]*rij[d]*r2inv; + } + } + } + m += 21; + } +------------------------------------------------------------------------- */ +} + + +/* ---------------------------------------------------------------------- + count angles and compute angle info on this proc + only count if 2nd atom is the one storing the angle + all atoms in interaction must be in group + all atoms in interaction must be known to proc + if bond is deleted or turned off (type <= 0) + do not count or count contribution + COMMENTED FOR NOW +---------------------------------------------------------------------- */ +void ComputeBorn::compute_angles() +{ +/* ---------------------------------------------------------------------- + int i,m,n,na,atom1,atom2,atom3,imol,iatom,atype,ivar; + tagint tagprev; + double delx1,dely1,delz1,delx2,dely2,delz2; + double rsq1,rsq2,r1,r2,cost; + double r1r2, r1r2inv; + double rsq1inv,rsq2inv,r1inv,r2inv,cinv; + double *ptr; + + double **x = atom->x; + tagint *tag = atom->tag; + int *num_angle = atom->num_angle; + tagint **angle_atom1 = atom->angle_atom1; + tagint **angle_atom2 = atom->angle_atom2; + tagint **angle_atom3 = atom->angle_atom3; + int **angle_type = atom->angle_type; + int *mask = atom->mask; + + int *molindex = atom->molindex; + int *molatom = atom->molatom; + Molecule **onemols = atom->avec->onemols; + + int nlocal = atom->nlocal; + int molecular = atom->molecular; + + // loop over all atoms and their angles + + Angle *angle = force->angle; + + int a,b,c,d,e,f; + double duang, du2ang; + double del1[3], del2[3]; + double dcos[6]; + double d2cos[21]; + double d2lncos[21]; + + // Initializing array for intermediate cos derivatives + // w regard to strain + for (i = 0; i < 6; i++) { + dcos[i] = 0; + } + for (i = 0; i < 21; i++) { + d2cos[i] = 0; + d2lncos[i] = 0; + } + + m = 0; + while (m < nvalues) { + for (atom2 = 0; atom2 < nlocal; atom2++) { + if (!(mask[atom2] & groupbit)) continue; + + if (molecular == 1) na = num_angle[atom2]; + else { + if (molindex[atom2] < 0) continue; + imol = molindex[atom2]; + iatom = molatom[atom2]; + na = onemols[imol]->num_angle[iatom]; + } + + for (i = 0; i < na; i++) { + if (molecular == 1) { + if (tag[atom2] != angle_atom2[atom2][i]) continue; + atype = angle_type[atom2][i]; + atom1 = atom->map(angle_atom1[atom2][i]); + atom3 = atom->map(angle_atom3[atom2][i]); + } else { + if (tag[atom2] != onemols[imol]->angle_atom2[atom2][i]) continue; + atype = onemols[imol]->angle_type[atom2][i]; + tagprev = tag[atom2] - iatom - 1; + atom1 = atom->map(onemols[imol]->angle_atom1[atom2][i]+tagprev); + atom3 = atom->map(onemols[imol]->angle_atom3[atom2][i]+tagprev); + } + + if (atom1 < 0 || !(mask[atom1] & groupbit)) continue; + if (atom3 < 0 || !(mask[atom3] & groupbit)) continue; + if (atype <= 0) continue; + + delx1 = x[atom1][0] - x[atom2][0]; + dely1 = x[atom1][1] - x[atom2][1]; + delz1 = x[atom1][2] - x[atom2][2]; + domain->minimum_image(delx1,dely1,delz1); + del1[0] = delx1; + del1[1] = dely1; + del1[2] = delz1; + + rsq1 = delx1*delx1 + dely1*dely1 + delz1*delz1; + rsq1inv = 1.0/rsq1; + r1 = sqrt(rsq1); + r1inv = 1.0/r1; + + delx2 = x[atom3][0] - x[atom2][0]; + dely2 = x[atom3][1] - x[atom2][1]; + delz2 = x[atom3][2] - x[atom2][2]; + domain->minimum_image(delx2,dely2,delz2); + del2[0] = delx2; + del2[1] = dely2; + del2[2] = delz2; + + rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2; + rsq2inv = 1.0/rsq2; + r2 = sqrt(rsq2); + r2inv = 1.0/r2; + + r1r2 = delx1*delx2 + dely1*dely2 + delz1*delz2; + r1r2inv = 1/r1r2; + // cost = cosine of angle + + cost = delx1*delx2 + dely1*dely2 + delz1*delz2; + cost /= r1*r2; + if (cost > 1.0) cost = 1.0; + if (cost < -1.0) cost = -1.0; + cinv = 1.0/cost; + + // The method must return derivative with regards + // to cos(theta)! + // Use the chain rule if needed: + // dU(t)/de = dt/dcos(t)*dU(t)/dt*dcos(t)/de + // with dt/dcos(t) = -1/sin(t) + angle->born(atype,atom1,atom2,atom3,duang,du2ang); + + // Voigt notation + // 1 = 11, 2 = 22, 3 = 33 + // 4 = 23, 5 = 13, 6 = 12 + a = 0; + b = 0; + c = 0; + d = 0; + for (i = 0; i<6; i++) { + a = albe[i][0]; + b = albe[i][1]; + dcos[i] = cost*(del1[a]*del2[b]+del1[b]*del2[a]*r1r2inv - + del1[a]*del1[b]*rsq1inv - del2[a]*del2[b]*rsq2inv); + } + for (i = 0; i<21; i++) { + a = albemunu[i][0]; + b = albemunu[i][1]; + c = albemunu[i][2]; + d = albemunu[i][3]; + e = albe[i][0]; + f = albe[i][1]; + d2lncos[i] = 2*(del1[a]*del1[b]*del1[c]*del1[d]*rsq1inv*rsq1inv + + del2[a]*del2[b]*del2[c]*del2[d]*rsq2inv*rsq2inv) - + (del1[a]*del2[b]+del1[b]*del2[a]) * + (del1[c]*del2[d]+del1[d]*del2[c]) * + r1r2inv*r1r2inv; + d2cos[i] = cost*d2lncos[i] + dcos[e]*dcos[f]*cinv; + values_local[m+i] += duang*d2cos[i] + du2ang*dcos[e]*dcos[f]; + } + } + } + m+=21; + } +------------------------------------------------------------------------- */ +} + +/* ---------------------------------------------------------------------- + count dihedrals on this proc + only count if 2nd atom is the one storing the dihedral + all atoms in interaction must be in group + all atoms in interaction must be known to proc + if flag is set, compute requested info about dihedral + COMMENTED FOR NOW +------------------------------------------------------------------------- */ + +void ComputeBorn::compute_dihedrals() +{ +/* ---------------------------------------------------------------------- + int i,m,n,nd,atom1,atom2,atom3,atom4,imol,iatom,dtype,ivar; + tagint tagprev; + double vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z,vb2xm,vb2ym,vb2zm; + double ax,ay,az,bx,by,bz,rasq,rbsq,rgsq,rg,ra2inv,rb2inv,rabinv; + double si,co,phi; + double *ptr; + + double **x = atom->x; + tagint *tag = atom->tag; + int *num_dihedral = atom->num_dihedral; + tagint **dihedral_atom1 = atom->dihedral_atom1; + tagint **dihedral_atom2 = atom->dihedral_atom2; + tagint **dihedral_atom3 = atom->dihedral_atom3; + tagint **dihedral_atom4 = atom->dihedral_atom4; + int *mask = atom->mask; + + int *molindex = atom->molindex; + int *molatom = atom->molatom; + Molecule **onemols = atom->avec->onemols; + + int nlocal = atom->nlocal; + int molecular = atom->molecular; + + // loop over all atoms and their dihedrals + + Dihedral *dihedral = force->dihedral; + + double dudih,du2dih; + int a,b,c,d,e,f; + double b1sq; + double b2sq; + double b3sq; + double b1b2; + double b1b3; + double b2b3; + double b1[3]; + double b2[3]; + double b3[3]; + + // Actually derivatives of the square of the + // vectors dot product. + double dmn[6]; + double dmm[6]; + double dnn[6]; + double d2mn[21]; + double d2mm[21]; + double d2nn[21]; + + double dcos[6]; + double d2cos[21]; + + for (i = 0; i < 6; i++) { + dmn[i] =0; + dmm[i] = 0; + dnn[i] = 0; + dcos[i] = 0; + } + for (i = 0; i < 21; i++) { + d2mn[i] = 0; + d2mm[i] = 0; + d2nn[i] = 0; + d2cos[i] = 0; + } + + m = 0; + while (m < nvalues) { + for (atom2 = 0; atom2 < nlocal; atom2++) { + if (!(mask[atom2] & groupbit)) continue; + + if (molecular == 1) nd = num_dihedral[atom2]; + else { + if (molindex[atom2] < 0) continue; + imol = molindex[atom2]; + iatom = molatom[atom2]; + nd = onemols[imol]->num_dihedral[iatom]; + } + + for (i = 0; i < nd; i++) { + if (molecular == 1) { + if (tag[atom2] != dihedral_atom2[atom2][i]) continue; + atom1 = atom->map(dihedral_atom1[atom2][i]); + atom3 = atom->map(dihedral_atom3[atom2][i]); + atom4 = atom->map(dihedral_atom4[atom2][i]); + } else { + if (tag[atom2] != onemols[imol]->dihedral_atom2[atom2][i]) continue; + tagprev = tag[atom2] - iatom - 1; + atom1 = atom->map(onemols[imol]->dihedral_atom1[atom2][i]+tagprev); + atom3 = atom->map(onemols[imol]->dihedral_atom3[atom2][i]+tagprev); + atom4 = atom->map(onemols[imol]->dihedral_atom4[atom2][i]+tagprev); + } + + if (atom1 < 0 || !(mask[atom1] & groupbit)) continue; + if (atom3 < 0 || !(mask[atom3] & groupbit)) continue; + if (atom4 < 0 || !(mask[atom4] & groupbit)) continue; + + // phi calculation from dihedral style harmonic + + // The method must return derivative with regards + // to cos(phi)! + // Use the chain rule if needed: + // dU(t)/de = dt/dcos(t)*dU(t)/dt*dcos(t)/de + // with dt/dcos(t) = -1/sin(t) + + dihedral->born(nd,atom1,atom2,atom3,atom4,dudih,du2dih); + + vb1x = x[atom1][0] - x[atom2][0]; + vb1y = x[atom1][1] - x[atom2][1]; + vb1z = x[atom1][2] - x[atom2][2]; + domain->minimum_image(vb1x,vb1y,vb1z); + b1[0] = vb1x; + b1[1] = vb1y; + b1[2] = vb1z; + b1sq = b1[0]*b1[0]+b1[1]*b1[1]+b1[2]*b1[2]; + + vb2x = x[atom3][0] - x[atom2][0]; + vb2y = x[atom3][1] - x[atom2][1]; + vb2z = x[atom3][2] - x[atom2][2]; + domain->minimum_image(vb2x,vb2y,vb2z); + b2[0] = vb2x; + b2[1] = vb2y; + b2[2] = vb2z; + b2sq = b2[0]*b2[0]+b2[1]*b2[1]+b2[2]*b2[2]; + + vb2xm = -vb2x; + vb2ym = -vb2y; + vb2zm = -vb2z; + domain->minimum_image(vb2xm,vb2ym,vb2zm); + + vb3x = x[atom4][0] - x[atom3][0]; + vb3y = x[atom4][1] - x[atom3][1]; + vb3z = x[atom4][2] - x[atom3][2]; + domain->minimum_image(vb3x,vb3y,vb3z); + b3[0] = vb3x; + b3[1] = vb3y; + b3[2] = vb3z; + b3sq = b3[0]*b3[0]+b3[1]*b3[1]+b3[2]*b3[2]; + + b1b2 = b1[0]*b2[0]+b1[1]*b2[1]+b1[2]*b2[2]; + b1b3 = b1[0]*b3[0]+b1[1]*b3[1]+b1[2]*b3[2]; + b2b3 = b2[0]*b3[0]+b2[1]*b3[1]+b2[2]*b3[2]; + + ax = vb1y*vb2zm - vb1z*vb2ym; + ay = vb1z*vb2xm - vb1x*vb2zm; + az = vb1x*vb2ym - vb1y*vb2xm; + bx = vb3y*vb2zm - vb3z*vb2ym; + by = vb3z*vb2xm - vb3x*vb2zm; + bz = vb3x*vb2ym - vb3y*vb2xm; + + rasq = ax*ax + ay*ay + az*az; + rbsq = bx*bx + by*by + bz*bz; + rgsq = vb2xm*vb2xm + vb2ym*vb2ym + vb2zm*vb2zm; + rg = sqrt(rgsq); + + ra2inv = rb2inv = 0.0; + if (rasq > 0) ra2inv = 1.0/rasq; + if (rbsq > 0) rb2inv = 1.0/rbsq; + rabinv = sqrt(ra2inv*rb2inv); + + co = (ax*bx + ay*by + az*bz)*rabinv; + si = rg*rabinv*(ax*vb3x + ay*vb3y + az*vb3z); + + if (co > 1.0) co = 1.0; + if (co < -1.0) co = -1.0; + phi = atan2(si,co); + + // above a and b are m and n vectors + // here they are integers indices + a = 0; + b = 0; + c = 0; + d = 0; + e = 0; + f = 0; + for (i = 0; i<6; i++) { + a = albe[i][0]; + b = albe[i][1]; + dmm[i] = 2*(b2sq*b1[a]*b1[b]+b1sq*b2[a]*b2[b] - b1b2*(b1[a]*b2[b]+b1[b]*b2[a])); + dnn[i] = 2*(b3sq*b2[a]*b2[b]+b2sq*b3[a]*b3[b] - b2b3*(b2[a]*b3[b]+b2[b]*b3[a])); + dmn[i] = b1b2*(b2[a]*b3[b]+b2[b]*b3[a]) + b2b3*(b1[a]*b2[b]+b1[b]*b2[a]) + - 2*(b1b3*b2[a]*b2[b]) - b2sq*(b1[a]*b3[b]+b1[b]*b3[a]); + dcos[i] = co*(rabinv*rabinv*dmn[i] - ra2inv*dmm[i] - rb2inv*dnn[i])/2.; + } + for (i = 0; i<21; i++) { + a = albemunu[i][0]; + b = albemunu[i][1]; + c = albemunu[i][2]; + d = albemunu[i][3]; + e = albe[i][0]; + f = albe[i][1]; + d2mm[i] = 4*(b1[a]*b1[b]*b2[c]*b2[d] + b1[c]*b1[d]*b2[a]*b2[b]) + - 8*(b1[a]*b2[b]+b1[b]*b2[a])*(b1[c]*b2[d]+b1[d]*b2[c]); + d2nn[i] = 4*(b2[a]*b2[b]*b3[c]*b3[d] + b2[c]*b2[d]*b3[a]*b3[b]) + - 8*(b2[a]*b3[b]+b2[b]*b3[a])*(b2[c]*b3[d]+b2[d]*b3[c]); + d2mn[i] = (b1[a]*b2[b]+b1[b]*b2[a])*(b2[c]*b3[d]+b2[d]*b3[c]) + + (b2[a]*b3[b]+b2[b]*b3[a])*(b1[c]*b2[d]+b1[d]*b2[d]) + - (b1[a]*b3[b]+b1[b]*b3[a])*(b2[c]*b2[d]+b2[c]*b2[d]) + - (b1[c]*b3[d]+b1[d]*b3[c])*(b2[a]*b2[b]+b2[a]*b2[b]); + d2cos[i] = co/2.*( + rabinv*rabinv*d2mn[i] + - rabinv*rabinv*rabinv*rabinv*dmn[e]*dmn[f] + + ra2inv*ra2inv*dmm[e]*dmm[f] + - ra2inv*d2mm[i] + + rb2inv*rb2inv*dnn[e]*dnn[f] + - rb2inv*d2nn[i] ); + values_local[m+i] += dudih*d2cos[i] + du2dih*dcos[e]*dcos[f]; + } + } + } + m+=21; + } +------------------------------------------------------------------------- */ +} diff --git a/src/compute_born.h b/src/compute_born.h new file mode 100644 index 0000000000..390922364c --- /dev/null +++ b/src/compute_born.h @@ -0,0 +1,125 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. + ------------------------------------------------------------------------- */ + +/*------------------------------------------------------------------------ + Contributing Authors : Germain Clavier (UCA) + --------------------------------------------------------------------------*/ + +#ifdef COMPUTE_CLASS + +ComputeStyle(born,ComputeBorn) + +#else + +#ifndef LMP_COMPUTE_BORN_H +#define LMP_COMPUTE_BORN_H + +#include "compute.h" + +namespace LAMMPS_NS { + + class ComputeBorn : public Compute { + public: + ComputeBorn(class LAMMPS *, int, char **); + virtual ~ComputeBorn(); + void init(); + void init_list(int, class NeighList *); + void compute_vector(); + + private: + + void compute_pairs(); + void compute_bonds(); + void compute_angles(); + void compute_dihedrals(); + void compute_impropers(); + + int me,nvalues; + int *which; + + int pairflag, bondflag, angleflag; + int dihedflag, impflag; + + int const albe[21][2] = { + {0,0}, // C11 + {1,1}, // C22 + {2,2}, // C33 + {1,2}, // C44 + {0,2}, // C55 + {0,1}, // C66 + {0,1}, // C12 + {0,2}, // C13 + {0,3}, // C14 + {0,4}, // C15 + {0,5}, // C16 + {1,2}, // C23 + {1,3}, // C24 + {1,4}, // C25 + {1,5}, // C26 + {2,3}, // C34 + {2,4}, // C35 + {2,5}, // C36 + {3,4}, // C45 + {3,5}, // C46 + {4,5} // C56 + }; + + int const albemunu[21][4] = { + {0,0,0,0}, // C11 + {1,1,1,1}, // C22 + {2,2,2,2}, // C33 + {1,2,1,2}, // C44 + {0,2,0,2}, // C55 + {0,1,0,1}, // C66 + {0,0,1,1}, // C12 + {0,0,2,2}, // C13 + {0,0,1,2}, // C14 + {0,0,0,2}, // C15 + {0,0,0,1}, // C16 + {1,1,2,2}, // C23 + {1,1,1,2}, // C24 + {1,1,0,2}, // C25 + {1,1,0,1}, // C26 + {2,2,1,2}, // C34 + {2,2,0,2}, // C35 + {2,2,0,1}, // C36 + {1,2,0,2}, // C45 + {1,2,0,1}, // C46 + {0,1,0,2} // C56 + }; + + double *values_local,*values_global; + double pos,pos1,dt,nktv2p,ftm2v; + class NeighList *list; + + }; + +} + +#endif +#endif + +/* ERROR/WARNING messages: + + E: Illegal ... command + + Self-explanatory. Check the input script syntax and compare to the + documentation for the command. You can use -echo screen as a + command-line option when running LAMMPS to see the offending line. + + E: ... style does not support compute born + + Some component of the force field (pair, bond, angle...) does not provide + a function to return the Born term contribution. + */ + diff --git a/src/dihedral.cpp b/src/dihedral.cpp index a0b6746a2c..201ec4151e 100644 --- a/src/dihedral.cpp +++ b/src/dihedral.cpp @@ -36,6 +36,7 @@ Dihedral::Dihedral(LAMMPS *lmp) : Pointers(lmp) allocated = 0; suffix_flag = Suffix::NONE; + born_enable = 0; maxeatom = maxvatom = maxcvatom = 0; eatom = nullptr; diff --git a/src/dihedral.h b/src/dihedral.h index 81daeef201..83c0a3b239 100644 --- a/src/dihedral.h +++ b/src/dihedral.h @@ -26,6 +26,7 @@ class Dihedral : protected Pointers { int allocated; int *setflag; int writedata; // 1 if writes coeffs to data file + int born_enable; double energy; // accumulated energy double virial[6]; // accumulated virial: xx,yy,zz,xy,xz,yz double *eatom, **vatom; // accumulated per-atom energy/virial @@ -55,6 +56,11 @@ class Dihedral : protected Pointers { virtual void read_restart_settings(FILE *){}; virtual void write_data(FILE *) {} virtual double memory_usage(); + virtual void born(int/*dtype*/, int/*at1*/, int/*at2*/, int/*at3*/, int /*at4*/, double& du, double& du2) + { + du = 0.0; + du2 = 0.0; + } protected: int suffix_flag; // suffix compatibility flag diff --git a/src/improper.cpp b/src/improper.cpp index 00033880f0..017388ddd7 100644 --- a/src/improper.cpp +++ b/src/improper.cpp @@ -34,6 +34,7 @@ Improper::Improper(LAMMPS *lmp) : Pointers(lmp) allocated = 0; suffix_flag = Suffix::NONE; + born_enable = 0; maxeatom = maxvatom = maxcvatom = 0; eatom = nullptr; diff --git a/src/improper.h b/src/improper.h index 1f88204a6b..94fcf62ab5 100644 --- a/src/improper.h +++ b/src/improper.h @@ -26,6 +26,7 @@ class Improper : protected Pointers { int allocated; int *setflag; int writedata; // 1 if writes coeffs to data file + int born_enable; double energy; // accumulated energies double virial[6]; // accumulated virial: xx,yy,zz,xy,xz,yz double *eatom, **vatom; // accumulated per-atom energy/virial @@ -55,6 +56,11 @@ class Improper : protected Pointers { virtual void read_restart_settings(FILE *){}; virtual void write_data(FILE *) {} virtual double memory_usage(); + virtual void born(int/*dtype*/, int/*at1*/, int/*at2*/, int/*at3*/, int /*at4*/, double& du, double& du2) + { + du = 0.0; + du2 = 0.0; + } protected: int suffix_flag; // suffix compatibility flag diff --git a/src/pair.cpp b/src/pair.cpp index f88c4e0972..2da47becca 100644 --- a/src/pair.cpp +++ b/src/pair.cpp @@ -58,6 +58,7 @@ Pair::Pair(LAMMPS *lmp) : Pointers(lmp) comm_forward = comm_reverse = comm_reverse_off = 0; single_enable = 1; + born_enable = 0; single_hessian_enable = 0; restartinfo = 1; respa_enable = 0; diff --git a/src/pair.h b/src/pair.h index 00e6734773..e793f588ae 100644 --- a/src/pair.h +++ b/src/pair.h @@ -51,6 +51,7 @@ class Pair : protected Pointers { int comm_reverse_off; // size of reverse comm even if newton off int single_enable; // 1 if single() routine exists + int born_enable; // 1 if born() routine exists int single_hessian_enable; // 1 if single_hessian() routine exists int restartinfo; // 1 if pair style writes restart info int respa_enable; // 1 if inner/middle/outer rRESPA routines @@ -158,6 +159,13 @@ class Pair : protected Pointers { return 0.0; } + virtual void born(int /*i*/, int /*j*/, int /*itype*/, int /*jtype*/, double /*rsq*/, + double /*factor_coul*/, double /*factor_lj*/, double& du, double& du2) + { + du = 0.0; + du2 = 0.0; + } + void hessian_twobody(double fforce, double dfac, double delr[3], double phiTensor[6]); virtual double single_hessian(int, int, int, int, double, double[3], double, double, diff --git a/src/pair_lj_cut.cpp b/src/pair_lj_cut.cpp index a9d45b9007..40d9611a02 100644 --- a/src/pair_lj_cut.cpp +++ b/src/pair_lj_cut.cpp @@ -40,6 +40,7 @@ using namespace MathConst; PairLJCut::PairLJCut(LAMMPS *lmp) : Pair(lmp) { respa_enable = 1; + born_enable = 1; writedata = 1; } @@ -680,6 +681,28 @@ double PairLJCut::single(int /*i*/, int /*j*/, int itype, int jtype, double rsq, /* ---------------------------------------------------------------------- */ +void PairLJCut::born(int /*i*/, int /*j*/, int itype, int jtype, double rsq, + double /*factor_coul*/, double factor_lj, + double &dupair, double &du2pair) +{ + double rinv,r2inv,r6inv,du,du2; + + r2inv = 1.0/rsq; + rinv = sqrt(r2inv); + r6inv = r2inv*r2inv*r2inv; + + // Reminder: lj1 = 48*e*s^12, lj2 = 24*e*s^6 + // so dupair = -forcelj/r = -fforce*r (forcelj from single method) + + du = r6inv * rinv * (lj2[itype][jtype] - lj1[itype][jtype]*r6inv); + du2 = r6inv * r2inv * (13*lj1[itype][jtype]*r6inv - 7*lj2[itype][jtype]); + + dupair = factor_lj*du; + du2pair = factor_lj*du2; +} + +/* ---------------------------------------------------------------------- */ + void *PairLJCut::extract(const char *str, int &dim) { dim = 2; diff --git a/src/pair_lj_cut.h b/src/pair_lj_cut.h index 15b1d92175..f94f4d98ee 100644 --- a/src/pair_lj_cut.h +++ b/src/pair_lj_cut.h @@ -40,6 +40,7 @@ class PairLJCut : public Pair { void write_data(FILE *); void write_data_all(FILE *); double single(int, int, int, int, double, double, double, double &); + void born(int, int, int, int, double, double, double, double &, double &); void *extract(const char *, int &); void compute_inner();