diff --git a/src/angle_harmonic.cpp b/src/angle_harmonic.cpp deleted file mode 100644 index d28afd76d6..0000000000 --- a/src/angle_harmonic.cpp +++ /dev/null @@ -1,264 +0,0 @@ -/* ---------------------------------------------------------------------- - LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - http://lammps.sandia.gov, Sandia National Laboratories - Steve Plimpton, sjplimp@sandia.gov - - Copyright (2003) Sandia Corporation. Under the terms of Contract - DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains - certain rights in this software. This software is distributed under - the GNU General Public License. - - See the README file in the top-level LAMMPS directory. -------------------------------------------------------------------------- */ - -#include -#include -#include "angle_harmonic.h" -#include "atom.h" -#include "neighbor.h" -#include "domain.h" -#include "comm.h" -#include "force.h" -#include "math_const.h" -#include "memory.h" -#include "error.h" - -using namespace LAMMPS_NS; -using namespace MathConst; - -#define SMALL 0.001 - -/* ---------------------------------------------------------------------- */ - -AngleHarmonic::AngleHarmonic(LAMMPS *lmp) : Angle(lmp) -{ - k = NULL; - theta0 = NULL; -} - -/* ---------------------------------------------------------------------- */ - -AngleHarmonic::~AngleHarmonic() -{ - if (allocated && !copymode) { - memory->destroy(setflag); - memory->destroy(k); - memory->destroy(theta0); - } -} - -/* ---------------------------------------------------------------------- */ - -void AngleHarmonic::compute(int eflag, int vflag) -{ - int i1,i2,i3,n,type; - double delx1,dely1,delz1,delx2,dely2,delz2; - double eangle,f1[3],f3[3]; - double dtheta,tk; - double rsq1,rsq2,r1,r2,c,s,a,a11,a12,a22; - - eangle = 0.0; - if (eflag || vflag) ev_setup(eflag,vflag); - else evflag = 0; - - double **x = atom->x; - double **f = atom->f; - int **anglelist = neighbor->anglelist; - int nanglelist = neighbor->nanglelist; - int nlocal = atom->nlocal; - int newton_bond = force->newton_bond; - - for (n = 0; n < nanglelist; n++) { - i1 = anglelist[n][0]; - i2 = anglelist[n][1]; - i3 = anglelist[n][2]; - type = anglelist[n][3]; - - // 1st bond - - delx1 = x[i1][0] - x[i2][0]; - dely1 = x[i1][1] - x[i2][1]; - delz1 = x[i1][2] - x[i2][2]; - - rsq1 = delx1*delx1 + dely1*dely1 + delz1*delz1; - r1 = sqrt(rsq1); - - // 2nd bond - - delx2 = x[i3][0] - x[i2][0]; - dely2 = x[i3][1] - x[i2][1]; - delz2 = x[i3][2] - x[i2][2]; - - rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2; - r2 = sqrt(rsq2); - - // angle (cos and sin) - - c = delx1*delx2 + dely1*dely2 + delz1*delz2; - c /= r1*r2; - - if (c > 1.0) c = 1.0; - if (c < -1.0) c = -1.0; - - s = sqrt(1.0 - c*c); - if (s < SMALL) s = SMALL; - s = 1.0/s; - - // force & energy - - dtheta = acos(c) - theta0[type]; - tk = k[type] * dtheta; - - if (eflag) eangle = tk*dtheta; - - a = -2.0 * tk * s; - a11 = a*c / rsq1; - a12 = -a / (r1*r2); - a22 = a*c / rsq2; - - f1[0] = a11*delx1 + a12*delx2; - f1[1] = a11*dely1 + a12*dely2; - f1[2] = a11*delz1 + a12*delz2; - f3[0] = a22*delx2 + a12*delx1; - f3[1] = a22*dely2 + a12*dely1; - f3[2] = a22*delz2 + a12*delz1; - - // apply force to each of 3 atoms - - if (newton_bond || i1 < nlocal) { - f[i1][0] += f1[0]; - f[i1][1] += f1[1]; - f[i1][2] += f1[2]; - } - - if (newton_bond || i2 < nlocal) { - f[i2][0] -= f1[0] + f3[0]; - f[i2][1] -= f1[1] + f3[1]; - f[i2][2] -= f1[2] + f3[2]; - } - - if (newton_bond || i3 < nlocal) { - f[i3][0] += f3[0]; - f[i3][1] += f3[1]; - f[i3][2] += f3[2]; - } - - if (evflag) ev_tally(i1,i2,i3,nlocal,newton_bond,eangle,f1,f3, - delx1,dely1,delz1,delx2,dely2,delz2); - } -} - -/* ---------------------------------------------------------------------- */ - -void AngleHarmonic::allocate() -{ - allocated = 1; - int n = atom->nangletypes; - - memory->create(k,n+1,"angle:k"); - memory->create(theta0,n+1,"angle:theta0"); - - memory->create(setflag,n+1,"angle:setflag"); - for (int i = 1; i <= n; i++) setflag[i] = 0; -} - -/* ---------------------------------------------------------------------- - set coeffs for one or more types -------------------------------------------------------------------------- */ - -void AngleHarmonic::coeff(int narg, char **arg) -{ - if (narg != 3) error->all(FLERR,"Incorrect args for angle coefficients"); - if (!allocated) allocate(); - - int ilo,ihi; - force->bounds(FLERR,arg[0],atom->nangletypes,ilo,ihi); - - double k_one = force->numeric(FLERR,arg[1]); - double theta0_one = force->numeric(FLERR,arg[2]); - - // convert theta0 from degrees to radians - - int count = 0; - for (int i = ilo; i <= ihi; i++) { - k[i] = k_one; - theta0[i] = theta0_one/180.0 * MY_PI; - setflag[i] = 1; - count++; - } - - if (count == 0) error->all(FLERR,"Incorrect args for angle coefficients"); -} - -/* ---------------------------------------------------------------------- */ - -double AngleHarmonic::equilibrium_angle(int i) -{ - return theta0[i]; -} - -/* ---------------------------------------------------------------------- - proc 0 writes out coeffs to restart file -------------------------------------------------------------------------- */ - -void AngleHarmonic::write_restart(FILE *fp) -{ - fwrite(&k[1],sizeof(double),atom->nangletypes,fp); - fwrite(&theta0[1],sizeof(double),atom->nangletypes,fp); -} - -/* ---------------------------------------------------------------------- - proc 0 reads coeffs from restart file, bcasts them -------------------------------------------------------------------------- */ - -void AngleHarmonic::read_restart(FILE *fp) -{ - allocate(); - - if (comm->me == 0) { - fread(&k[1],sizeof(double),atom->nangletypes,fp); - fread(&theta0[1],sizeof(double),atom->nangletypes,fp); - } - MPI_Bcast(&k[1],atom->nangletypes,MPI_DOUBLE,0,world); - MPI_Bcast(&theta0[1],atom->nangletypes,MPI_DOUBLE,0,world); - - for (int i = 1; i <= atom->nangletypes; i++) setflag[i] = 1; -} - -/* ---------------------------------------------------------------------- - proc 0 writes to data file -------------------------------------------------------------------------- */ - -void AngleHarmonic::write_data(FILE *fp) -{ - for (int i = 1; i <= atom->nangletypes; i++) - fprintf(fp,"%d %g %g\n",i,k[i],theta0[i]/MY_PI*180.0); -} - -/* ---------------------------------------------------------------------- */ - -double AngleHarmonic::single(int type, int i1, int i2, int i3) -{ - double **x = atom->x; - - double delx1 = x[i1][0] - x[i2][0]; - double dely1 = x[i1][1] - x[i2][1]; - double delz1 = x[i1][2] - x[i2][2]; - domain->minimum_image(delx1,dely1,delz1); - double r1 = sqrt(delx1*delx1 + dely1*dely1 + delz1*delz1); - - double delx2 = x[i3][0] - x[i2][0]; - double dely2 = x[i3][1] - x[i2][1]; - double delz2 = x[i3][2] - x[i2][2]; - domain->minimum_image(delx2,dely2,delz2); - double r2 = sqrt(delx2*delx2 + dely2*dely2 + delz2*delz2); - - double c = delx1*delx2 + dely1*dely2 + delz1*delz2; - c /= r1*r2; - if (c > 1.0) c = 1.0; - if (c < -1.0) c = -1.0; - - double dtheta = acos(c) - theta0[type]; - double tk = k[type] * dtheta; - return tk*dtheta; -} diff --git a/src/bond_harmonic.cpp b/src/bond_harmonic.cpp deleted file mode 100644 index f795610b37..0000000000 --- a/src/bond_harmonic.cpp +++ /dev/null @@ -1,215 +0,0 @@ -/* ---------------------------------------------------------------------- - LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - http://lammps.sandia.gov, Sandia National Laboratories - Steve Plimpton, sjplimp@sandia.gov - - Copyright (2003) Sandia Corporation. Under the terms of Contract - DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains - certain rights in this software. This software is distributed under - the GNU General Public License. - - See the README file in the top-level LAMMPS directory. -------------------------------------------------------------------------- */ - -#include -#include -#include -#include "bond_harmonic.h" -#include "atom.h" -#include "neighbor.h" -#include "domain.h" -#include "comm.h" -#include "force.h" -#include "memory.h" -#include "error.h" - -using namespace LAMMPS_NS; - -/* ---------------------------------------------------------------------- */ - -BondHarmonic::BondHarmonic(LAMMPS *lmp) : Bond(lmp) -{ - reinitflag = 1; -} - -/* ---------------------------------------------------------------------- */ - -BondHarmonic::~BondHarmonic() -{ - if (allocated && !copymode) { - memory->destroy(setflag); - memory->destroy(k); - memory->destroy(r0); - } -} - -/* ---------------------------------------------------------------------- */ - -void BondHarmonic::compute(int eflag, int vflag) -{ - int i1,i2,n,type; - double delx,dely,delz,ebond,fbond; - double rsq,r,dr,rk; - - ebond = 0.0; - if (eflag || vflag) ev_setup(eflag,vflag); - else evflag = 0; - - double **x = atom->x; - double **f = atom->f; - int **bondlist = neighbor->bondlist; - int nbondlist = neighbor->nbondlist; - int nlocal = atom->nlocal; - int newton_bond = force->newton_bond; - - for (n = 0; n < nbondlist; n++) { - i1 = bondlist[n][0]; - i2 = bondlist[n][1]; - type = bondlist[n][2]; - - delx = x[i1][0] - x[i2][0]; - dely = x[i1][1] - x[i2][1]; - delz = x[i1][2] - x[i2][2]; - - rsq = delx*delx + dely*dely + delz*delz; - r = sqrt(rsq); - dr = r - r0[type]; - rk = k[type] * dr; - - // force & energy - - if (r > 0.0) fbond = -2.0*rk/r; - else fbond = 0.0; - - if (eflag) ebond = rk*dr; - - // apply force to each of 2 atoms - - if (newton_bond || i1 < nlocal) { - f[i1][0] += delx*fbond; - f[i1][1] += dely*fbond; - f[i1][2] += delz*fbond; - } - - if (newton_bond || i2 < nlocal) { - f[i2][0] -= delx*fbond; - f[i2][1] -= dely*fbond; - f[i2][2] -= delz*fbond; - } - - if (evflag) ev_tally(i1,i2,nlocal,newton_bond,ebond,fbond,delx,dely,delz); - } -} - -/* ---------------------------------------------------------------------- */ - -void BondHarmonic::allocate() -{ - allocated = 1; - int n = atom->nbondtypes; - - memory->create(k,n+1,"bond:k"); - memory->create(r0,n+1,"bond:r0"); - - memory->create(setflag,n+1,"bond:setflag"); - for (int i = 1; i <= n; i++) setflag[i] = 0; -} - -/* ---------------------------------------------------------------------- - set coeffs for one or more types -------------------------------------------------------------------------- */ - -void BondHarmonic::coeff(int narg, char **arg) -{ - if (narg != 3) error->all(FLERR,"Incorrect args for bond coefficients"); - if (!allocated) allocate(); - - int ilo,ihi; - force->bounds(FLERR,arg[0],atom->nbondtypes,ilo,ihi); - - double k_one = force->numeric(FLERR,arg[1]); - double r0_one = force->numeric(FLERR,arg[2]); - - int count = 0; - for (int i = ilo; i <= ihi; i++) { - k[i] = k_one; - r0[i] = r0_one; - setflag[i] = 1; - count++; - } - - if (count == 0) error->all(FLERR,"Incorrect args for bond coefficients"); -} - -/* ---------------------------------------------------------------------- - return an equilbrium bond length -------------------------------------------------------------------------- */ - -double BondHarmonic::equilibrium_distance(int i) -{ - return r0[i]; -} - -/* ---------------------------------------------------------------------- - proc 0 writes out coeffs to restart file -------------------------------------------------------------------------- */ - -void BondHarmonic::write_restart(FILE *fp) -{ - fwrite(&k[1],sizeof(double),atom->nbondtypes,fp); - fwrite(&r0[1],sizeof(double),atom->nbondtypes,fp); -} - -/* ---------------------------------------------------------------------- - proc 0 reads coeffs from restart file, bcasts them -------------------------------------------------------------------------- */ - -void BondHarmonic::read_restart(FILE *fp) -{ - allocate(); - - if (comm->me == 0) { - fread(&k[1],sizeof(double),atom->nbondtypes,fp); - fread(&r0[1],sizeof(double),atom->nbondtypes,fp); - } - MPI_Bcast(&k[1],atom->nbondtypes,MPI_DOUBLE,0,world); - MPI_Bcast(&r0[1],atom->nbondtypes,MPI_DOUBLE,0,world); - - for (int i = 1; i <= atom->nbondtypes; i++) setflag[i] = 1; -} - -/* ---------------------------------------------------------------------- - proc 0 writes to data file -------------------------------------------------------------------------- */ - -void BondHarmonic::write_data(FILE *fp) -{ - for (int i = 1; i <= atom->nbondtypes; i++) - fprintf(fp,"%d %g %g\n",i,k[i],r0[i]); -} - -/* ---------------------------------------------------------------------- */ - -double BondHarmonic::single(int type, double rsq, int /*i*/, int /*j*/, - double &fforce) -{ - double r = sqrt(rsq); - double dr = r - r0[type]; - double rk = k[type] * dr; - fforce = 0; - if (r > 0.0) fforce = -2.0*rk/r; - return rk*dr; -} - -/* ---------------------------------------------------------------------- - Return ptr to internal members upon request. ------------------------------------------------------------------------- */ -void *BondHarmonic::extract( char *str, int &dim ) -{ - dim = 1; - if( strcmp(str,"kappa")==0) return (void*) k; - if( strcmp(str,"r0")==0) return (void*) r0; - return NULL; -} - -